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Foreword 



Commodity prices are an important input factor for capital investment decisions, 
for instance when it comes to valuation of oil fields, mines or chemical plant 
configurations to name but a few areas of application. One investment appraisal 
method which allows to explicitly quantify managerial flexibility is real options 
valuation which requires the proper specification of the underlying - typically the 
commodity price - in terms of a stochastic process. Most of the literature relies on 
the assumption that commodity prices are characterised by Geometric Brownian 
Motion (GBM) which entails mathematical tractability. Max Schone, as a starting 
point for his Master’s Thesis, questions this premise and sets out to explore alter- 
native stochastic processes to model commodity price processes providing a better 
fit to historical data. 

Schone devises in his thesis a calibration scheme which allows him to assess the 
fit of the investigated processes to a basket of commodities. He finds that the alter- 
natives to GBM provide a better fit and singles out the variance gamma process as 
the best compromise between goodness of fit and model complexity. He then goes 
on to demonstrate that the more appropriate specification of the commodity price 
process for the purpose of real options valuation may affect both project values by 
as much as 20% and the optimal time to invest. 

This thesis therefore serves for the practitioner in capital investment valuation as 
an excellent source for model selection and for the researcher as a solid starting 
point for enhanced commodity price modelling. 



Prof. Dr. Stefan Spinier 
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1 Introduction 



The success or failure of businesses is determined by the ability of management 
to systematically seize value creating investment opportunities. To this end, an 
accurate investment valuation is of paramount importance as it resembles the only 
broad and long-term oriented criterion for investment decision making (Brealey, 
Allen, & Myers, 2011; Roller, Wessels, & Goedhart, 2010). 

While Discounted Cash Flow (DCF) and Net Present Value (NPV) analyses are 
the easiest and most widespread methods to obtain an estimate of project value, 
they fail to account for the value of managerial flexibility inherent in most invest- 
ment decisions (Trigeorgis, 1996). Accordingly, Copeland and Antikarov (2003) 
argue that ’'NPV systematically undervalues every project” (p. 5), and thereby re- 
inforce the findings of Bhappu and Guzman (1995) and K. Roberts and Weitzman 
(as cited in Schulmerich, 2010), who contend that even negative NPV projects 
may be worthwhile to undertake in cases of sequential decision making. Since, 
very few managerial decisions are completely irreversible or alterable (Copeland 
& Antikarov, 2003), the idea of Dixit and Pindyck (1994) to capture the arising 
value of flexibility using options theory, revolutionised the world of corporate fi- 
nance and laid the foundation for the concept of real option valuation. 

Although, this option-based framework unites indisputable benefits in capturing 
the full ramifications of uncertainty and flexibility in investment valuation, it comes 
at the cost of higher complexity. In particular, to capture the value arising from 
management’s ability to revise decisions contingent on the future state of the world, 
it is of central importance to make realistic assumptions on the probability distri- 
bution of future events. If this was not the case, one would account for the value 
of contingent actions based on outcomes that are unrealistic in magnitude or fre- 
quency in the first place. Clearly, under such a scenario valuation results cannot be 
meaningful and reliable. 

In real options analysis, the role of simulating future events accrues to a stochastic 
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process that drives the evolution of prices, costs, or output (Copeland & Antikarov, 
2003). Moreover, one or more of these quantities are often driven by uncertain fu- 
ture commodity prices as these represent the selling price for natural resource pro- 
ducers such as mining companies (Brennan & Schwartz, 1985), petroleum firms 
(Siegel, Smith, & Paddock, 1987), or agricultural enterprises (Nissanke, 2012) and 
in other cases may account for a significant share of production costs at indus- 
trial firms. While not only the investment volume in those sectors is high, 1 it has 
also become increasingly challenging for businesses to cope with the uncertainty 
of volatile commodity prices in managerial decision making (Schuh, Strohmer, & 
Zetterqvist, 201 1). As a result, the concept of option based valuation is as timely as 
ever but, at the same time, relies critically on the ability of managers and investors 
to model the stochastic evolution of future commodity prices appropriately. How- 
ever, as outlined in more detail in the following subsection ’’empirical studies on 
commodity markets are sparse” (Brooks & Prokopczuk, 2013, p. 528) and often 
provide little practical guidance for the choice of stochastic models in commodity- 
linked capital investment valuation. 

1 . 1 Problem and objective 

The choice of stochastic process matters in real option valuation (Laughton & Ja- 
coby, 1993; Tsekrekos, Shackleton, & Wojakowski, 2012), however, the spectrum 
of stochastic processes that has been considered in this context is narrow com- 
pared to the array of research available on financial options (Brooks & Prokopczuk, 
2013). In fact, the majority of studies on option-based capital investment valua- 
tion relies on two types of processes (Ozorio, Bastian-Pinto, & Brandao, 2011). 
First, Geometric Brownian Motion (GBM), as suggested by Black and Scholes 
(1973) and Cox, Ross, and Rubinstein (1979), is used, among others, by Brennan 
and Schwartz (1985), Siegel et al. (1987), and Copeland and Antikarov (2003). 



According to Esty (2004), project financing in the natural resource industry alone accounted for more than $34 
billion in 2001. 
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Second, variations of mean-reverting models, are applied, for instance, by Gibson 
and Schwartz (1990), Ross (1997), Schwartz (1997), Schwartz and Smith (2000), 
Manoliu and Tompaidis (2002), and Casassus and Collin-Dufresne (2005). 2 De- 
spite thorough justification in each case and general acceptance of several of these 
processes in the literature, there are likewise a number reasons to question their 
applicability to capital investment valuation. In particular, the assumptions of nor- 
mally distributed spot price returns has led to criticism of GBM in the context 
of both, financial options (Brooks & Prokopczuk, 2013; Cont, 2001; Schoutens, 
2003) and real options (Lund, 1992; Metcalf & Hassett, 1995) and also the idea 
of mean reversion in commodity prices appears debatable. While it is consistent 
with economic intuition 3 and empirical evidence from the term structure of futures 
prices (Bessembinder, Coughenour, Seguin, & Smoller, 1995; Casassus & Collin- 
Dufresne, 2005), Dixit and Pindyck (1994) find similar support in historical price 
series only for long time spans in excess of 40 years. Moreover, Pindyck (1999) 
fails to reject a unit root 4 in coal and natural gas data over the period from 1875 
to 1996 and requires more than 100 years of data to reach the conclusion of mean 
reversion in crude oil prices. As, however, investment horizons in capital budget- 
ing are often significantly shorter, 5 it is unclear if mean-reverting models resemble 
a better choice than GBM in the valuation of commodity-linked contingent claims 
(Pindyck, 1999). 

Besides this ambiguity surrounding the two mainly considered types of stochas- 
tic models in a capital investment context, no single comprehensive and practical 
study is available on the subject. While some of the previously mentioned research 

2 Note that in some models (e.g. the two-factor model of Gibson and Schwartz (1990)) mean reversion is induced 
through a positive correlation between the spot price and a convenience yield process. In other words, when spot prices 
rise, an increase in the convenience yield will in turn decrease the drift of the spot price process which ’’will have a 
mean reversion effect on the spot commodtiy price” (Schwartz & Miltersen, 1998, p. 34). This approach has to be 
distinguished from processes of the Omstein-Uhlenbeck type (e.g. the one-factor model of Schwartz (1997)). 

3 When commodity prices are high, supply should increase as higher cost producers begin to enter the market and 
thereby drive prices back down. Conversely, when commodity prices are low, higher cost producers will be driven out 
of the market leading to a decrease in supply and rising prices. Over time, prices are expected to revert to long-run 
total marginal cost Pindyck (1999). 

4 The concept of unit root tests is discussed in section 3.1. 

5 Siegel et al. (1987), for instance, report that the relinquishment of a real world offshore oil property lease is (only) 
ten years. 
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has focused on the development of commodity price models to obtain a better fit 
to market prices of futures contracts, 6 these reports were mainly grounded on gen- 
erally accepted theoretical arguments on mean reversion and convenience yield 
dynamics, 7 but did not pursue the econometric verification of mean reversion in 
past data. Those studies, however, that undertook an empirical investigation did 
not analyse the merits of alternative price processes that are more consistent with 
identified scarce empirical support for mean reversion. 8 Furthermore, the range of 
commodities covered in these analyses is often narrow so that general conclusions 
on the relative appropriateness of stochastic processes in different markets are dif- 
ficult to derive, given the heterogeneity of commodity price dynamics documented 
by Kat and Oomen (2007). 

What follows from the previously presented arguments can be summarised in three 
main objectives for this paper. First, a preferably broad empirical analysis of com- 
modity price dynamics is necessary to shed additional light on the question of re- 
turn distributions and mean reversion as it will allow me to verify the appropriate- 
ness of currently popular GBM and mean-reverting models in capital investment 
valuation. Second, potential alternative stochastic processes need to be consid- 
ered that are more consistent with empirical evidence. Third, it must be clarified 
whether the choice of different stochastic models has a significant impact on valu- 
ation results and optimal investment decision rules, so that we can gauge the extent 
to which model selection deserves careful attention in practice. 

1.2 Course of investigation 

The remainder of this thesis is organised as follows. Section two presents the 
dataset. Next, an empirical assessment of mean reversion in commodity prices and 
an investigation of the normality of return distributions is provided in part three. 



6 See for instance Schwartz (1997), Schwartz and Smith (2000), or Casassus and Collin-Dufresne (2005). 

7 See Brennan (1991) for further reference on the role of convenience yield modelling. 

8 See Dixit and Pindyck (1994) and Pindyck (1999) for a discussion on mean reversion in historical commodity 
prices. 
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Once we have a better understanding of empirical price dynamics, it is investigated 
what characteristics a stochastic process should possess in order to realistically re- 
produce empirical price dynamics. Accordingly, a range of potential models with 
different properties is outlined in the first half of section four. Subsequently, a 
practical and flexible calibration method is proposed that facilitates the assessment 
of relative goodness of fit of each process to different commodity data. Based 
on the calibrated price models it is analysed in section five to what extent pro- 
cess choice influences the valuation of a stylised capital investment project with 
multiple inherent options. Section six concludes. 
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2 Data 



The task of setting up an appropriate dataset for my purposes involves decisions 
on the commodities to include, whether to use short-term future contracts or spot 
price data, what time period to analyse, and what data frequency to use. 

Following Marshall, Nguyen, and Visaltanachoti (2012), I focus on the compo- 
nents of the S&P Goldman Sachs Commodity Index (S&P GSCI). As this index 
comprises 24 important commodities spanning energy, industrial metal, precious 
metal, agriculture, and livestock, it suits my desire to obtain conclusions of prefer- 
ably broad, cross-sectoral applicability and to account for the heterogeneity of 
different commodities (Kat & Oomen, 2007). 

With respect to the type of price quote to use, I would generally prefer to use spot 
instead of future prices for my purposes, however, Fama and French (1987) sug- 
gest that ’’good spot-price data are not available for most commodities” (p. 57) 
and Schwartz (1997) raises similar concerns about their liquidity. As a result, 
these authors propose the use of short-term futures data instead. Since this view 
is not shared by Brooks and Prokopczuk (2013) and would generally complicate 
the analysis and interpretation of results throughout the paper, I conjecture that in 
the wake of ’’tremendous growth in commodity [...] markets” (Tsekrekos et ah, 
2012, p. 543), these earlier findings are not necessarily reflective of spot price data 
quality in the more recent past. 9 Since the use of justifiable data can be pivotal for 
the results of any empirical study (Gujarati, 2003), I feel it is necessary to assess 
the relative data quality in spot and futures markets based on a short, back-of-the- 
envelope liquidity calculation. 

For this task, the simple but computationally feasible 10 Zeros measure turns out 
useful (Lesmond, Ogden, & Trzcinka, 1999). For a given underlying, this met- 

9 Fama and French (1987) and Schwartz (1997) use datasets for the periods 1966 - 1984 and 1985 - 1995, respec- 
tively. 

10 Marshall et al. (2012) advise the Amihud (2002) or Effective Tick measure developed by (Goyenko, Holden, & 
Trzcinka, 2009) to calculate the liquidity of a security, however, these measures rely on trading volume data or price 
spreads, which are not consistently available from Thomson Reuters for spot price data of S&P GSCI components. 

6 
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Table 1: 

Overview of commodities in the dataset 

BBL: Barrel, BSH: Bushel, FED: U.S. Federal Reserve, ICCO: International Cocoa Organization, ISO: In- 
ternational Sugar Organization, LBM: London Bullion Market, LME: London Metal Exchange, MT: Metric 
Ton, OZT: Troy Ounce, TR: Thomson Reuters, USDA: U.S. Department of Agriculture. Note that the data 
is presented following the logic of four distinctive areas: Energy, industrial metals, precious metals, and 
agriculture. The sample period for all items is from 02/08/1993 to 30/12/2013. All data are in US$. 



Name 


Code 


Full description 


Unit 


Source 


Crude oil 


CRUDOIL 


Crude Oil-WTI Spot Cushing 


BBL 


TR 


Aluminum 


LAHCASH 


LME-Aluminium 99.7% Cash 


MT 


LME 


Copper 


LCPCASH 


LME-Copper Grade A Cash 


” 


” 


Lead 


LEDCASH 


LME-Lead Cash 


” 


” 


Nickel 


LNICASH 


LME-Nickel Cash 


” 


” 


Zinc 


LZZCASH 


LME-SHG Zinc 99.995% Cash 


” 


” 


Gold 


GOLDBLN 


Gold Bullion LBM 


OZT 


LBM 


Silver 


SLVCASH 


Silver Fix LBM Cash 


” 


” 


Wheat 


WHEATHD 


Wheat, No.2 Hard (Kansas) 


BSH 


USDA 


Corn 


COTSCIL 


Corn US No.2 South Central IL 


” 


” 


Soybeans 


SOYADSC 


Yellow Soybn US NO.l Sth Dvprt 


” 


” 


Sugar 


WSUGDLY 


Raw Sugar-ISA 


LB 


ISO 


Coffee 


COFBRAZ 


Coffee-Brazilian (NY) 


” 


TR 


Cocoa 


COCINUS 


Cocoa-ICCO 


MT 


ICCO 



ric relies on the proven negative relationship between the number of trading days 
with zero returns and liquidity in a given time interval. Hence, we can gain a rudi- 
mentary understanding of the relative liquidity and data quality of spot and future 
prices by simply comparing the average number of zero-return trading days per in- 
terval in both markets. As directly comparable spot and 3-months future prices are 
only available from the London Metal Exchange (LME) for the industrial metals 
part of the S&P GSCI index, this short liquidity analysis is somewhat limited in 
scope, but, nonetheless, the exercise reveals that over the more recent past LME 
data quality is, in fact, higher for spot prices than future prices. By investigating 
the development of Zeros over time, this conclusion can be reconciled with previ- 
ous research as it shows unsatisfactory spot price data quality prior to July 1993, 
but a structural improvement thereafter. To conserve space, the graphical results of 
this analysis are provided in appendix Al. 

Following the previous argument, I use spot price data for all commodities with 
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satisfactory data availability for the period 02/08/1993 - 30/12/2013. With respect 
to the data frequency, I use weekly quotes in the majority of calculations, as I sus- 
pect these to be more reflective of decision making horizons in capital investment 
decisions. However, as the analysis of empirical time series properties and calibra- 
tion tasks can sometimes benefit from the information contained in more frequent 
daily returns, I follow Kat and Oomen (2007) and Brooks and Prokopczuk (2013) 
and resort to daily quotes in some cases. 

An overview of the 14 commodities in the dataset is given in table 1. A com- 
plementary graph showing the historical price evolution for the four commodity 
sectors and the S&P 500 equity index is given in appendix A2. A table with all 
non-trading dates that have been erased from the dataset for the previous and all 
subsequent analyses is given in appendix A3. Please note that it is due to limited 
availability of high quality spot price data for energy commodities and live stock 
that only 14 out of the 24 S&P GSCI components are used. 



3 Empirical analysis 



As a first step to a better understanding of commodity price dynamics, it is use- 
ful to exploit the rich set of information contained in historical price series via a 
range of econometric and statistical tests. Although, one may be sceptical that a 
backward looking analysis can yield valuable insights for the decision of how to 
model seemingly random price movements in the future, more than half a cen- 
tury of empirical studies have revealed that the statistical properties of such price 
variations are indeed common across numerous asset classes in different markets 
and time periods (Cont, 2001; Mantegna & Stanley, 2000). For my purposes, it 
is logical to split the empirical analysis into two main blocks. First, it is analysed 
whether commodity prices are mean-reverting. The answer to this question is a 
comer stone in this thesis, as it resembles a central assumption to distinguish be- 
tween stochastic processes. Second, we will assess whether commodity returns are 
normally distributed. 

3.1 Stationarity of prices 

Following from the previous discussion, the question of mean-reverting prices is 
the principal distinctive assumption between commonly used GBM and mean- 
reverting processes. Broadly speaking, a stochastic process is mean-reverting (in- 
terchangeably we may refer to such a process as stationary or one without unit 
root) if its mean and variance are constant over time and the covariance between 
two observations depends only on the lag between them but not on the time when 
it is computed (Gujarati, 2003). 11 On the contrary, a non-stationary price process 
is characterised by a growing and unlimited variance over time that allows prices 
to rise without bound. To see more explicitly how the property of mean reversion 



11 In the literature, this form of stationarity is known as weak stationarity. In the case of a strictly stationary series 
all moments of the probability distribution are invariant over time (Cryer & Chan, 2008). 
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enters a stochastic process let us consider 



S t = QS t . l +e„ (3.1.1) 

where S t denotes, for instance, the price of a commodity at time t, 0 the autore- 
gressive coefficient, and 8, a white noise process 12 . In this framework, the absolute 
value of 0 is governing the mechanics of mean reversion. In particular, we distin- 
guish between two cases. First, the non-stationary Random Walk (RW) process, 
where |0| = 1 and, second, the stationary AR(1) process where |0| < l. 13 As a 
deeper understanding of the link between 0, process variance, and stationarity is 
very beneficial for subsequent statistical testing, let us study these cases in greater 
detail. To begin with, it is useful to apply repeated backward substitution to rewrite 
our expression for .S', as 

St = 0 [0 [eSf — 3 +£f— 2 ] +E ( _i] + 8 , = 0 V 3 + 9“£i — 2 + 0£f-l +£;• (3.1.2) 

Considering now the case where 0 = 1, it can be seen that if we continue to substi- 
tute S and set the initial price to zero, the above equation simplifies to 

t 

S/ = £e«. (3-1.3) 

;=o 

This illustrates that under a RW model, the price of e.g. a commodity at an arbitrary 
point in time t is given by a sequence of random shocks, where the impact of every 
individual shock never fades away over time. Intuitively, this also explains why the 
range of possible outcomes and, hence, the variance of this non-stationary process 
should grow over time. To see more formally why this is the case let use consider 
the Mean Squared Error (MSE) of a price forecast in a RW process .v periods ahead. 



12 White noise refers to a series of serially uncorrelated innovations with a mean of zero and a constant non-negative 
variance (Gujarati, 2003). 

13 Note that when 0 = 0, (3.1.1) collapses to a (stationary) white noise process. 
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Hence, we are interested in the quantity 



MS£. s = £[(S /+ ,-S f+i ) 2 ], (3.1.4) 

where E [S I+J ] is the conditional expectation of S t s periods ahead and E [S r+J ] the 
unconditional expectation. Since S t can be expressed according to (3.1.3), we can 
write for the case with an initial price different from zero 



MSE s = E 



(St + ^ E t+i ~ St+s) 2 ■ 

i=l 



However, since the best unconditional estimate for the future price in this Marko- 
vian setting S /+5 is simply today’s price S t and var (e f ) = a 2 , we get 



MSE s = E 




(3.1.5) 



Hence, the variance of this non-stationary process grows linearly in time (s) and 
the volatility at the square root of time. 

In the second case where |0| < 1, we can observe in (3.1.2) that the influence of 
the past on current prices is declining in time, thereby restraining the process to 
fluctuate around its unconditional mean (zero in our example). To quantify this 
dynamic we are, again, interested in an analytic expression for the variance as a 
function of time. To compute MSE s recap that cov (e t , Er+j) = 0 (for s ^ 0) and that 
the unconditional expectation of an AR( 1 ) process is given by E [S /+i | 5) .... . So] = 
0'S,. As a first step, let us consider 



MS£ 3 =£[(S /+ 3-S, +3 ) 2 ] 




( Q 2 S t + 0“£r+l + 0£f+2 + £/+3 
(l+0 2 + 0 4 ) , 




2 



where the conditional expectation S, +3 is obtained by adapting the time index of 
(3.1.2). Computing MSE also at some other values of s, one will soon recognise 
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the pattern that leads to the following generalisation for s periods: 

MSE s = a 2 (l + 0 2 + • • • + 0 2i “ 2 ) . 

Applying now a simple mathematical trick, we multiply the previous expression 
by 9 2 to get 

Q 2 MSE s = a 2 (e 2 + • • • + 0 2 ^ 

and subtract the result again from the previous equation. Simplifying yields 

MSE s ( 1 - 0 2 ) = 0 2 ( 1 + 0 2s ) 

7 ( 1 + 0 2s ) 

MSE s = g2 ^ 1 _ 02 ■ (3-1-6) 

To see that this conditional variance converges to the unconditional variance as an 
upper bound as s increases, we can take the limit of (3.1.6) and obtain 

lim =y[S,] =T ?l (3.1,7) 

While it took us some steps to work through this exercise, it will soon become ap- 
parent how the role of process variance and the coefficient 0 are intimately related 
to subsequent stationarity tests and the underlying principles of common stochas- 
tic processes. In fact, all stochastic diffusion models are closely related to (3.1.1). 
GBM, for instance, is roughly speaking a RW with drift in continuous time that 
restricts prices to remain positive over time (Bjork, 2009). On the other hand, 
Ornstein-Uhlenbeck processes such as the popular Vasicek (1977) model are, in 
some variation, the continuous time limit of the simple AR(1) process. To round 
off the introductory discussion of stationarity, please consider fig. 1 for a graphical 
illustration of the impact 0 has on simulated price trajectories. Turning now to the 
actual question of this subsection, namely whether historical commodity prices are 
mean-reverting, we will first directly test the value of 0 and, second, analyse the 
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A: Random walk (0= 1) B: AR(1) process (0= 0.8) 





Figure 1: Implications of stationarity. Panel A illustrates the non-stationary behaviour of a random walk 
without drift and standard deviation a = 0.15. Panel B shows a stationary AR(1) process with unconditional 
long-term mean of zero and an autoregressive coefficient 0 = 0.8 as shown above. Both processes are 
computed with the identical pseudo random numbers over 100 simulation steps. 



development of empirical price variance over time to make inferences on station- 
arity. 

A common test of the hypothesis Hq : |0| = 1 against the alternative H\ : |0| < 1 is 
the Augmented Dickey-Fuller (ADF) test developed by Dickey and Fuller (1979, 
1981). The test involves a regression of the following kind: 

St = a + 8r + 9S)_i + PiAS)_i + P2AS)_2 4 F P pAS t - p 4- £f , (3.1.8) 

where S t denotes the natural logarithm of the commodity price. 14 PAS are added 
to control for serial correlation among the innovations (e) and AS = S t — S t - 1 . The 
lag order is given by p and a and 5 1 denote a drift and deterministic time trend, 
respectively. While drift and trend can be excluded, they allow to test the unit root 
hypothesis more accurately for different data characteristics. In other words, if the 
data series shows a linear trend, the term 5 1 can be added to control for the trend 
and test the hypothesis of trend stationarity. If the data instead exhibits no linear 
trend but seems to have a non-zero mean, the drift term a should be added. 15 
To account for the possibility of these data characteristics, I run the regression 

14 Taking logs is often advised in the literature if a series shows exponential growth to establish a better fit in linear 
regressions. See, for instance, Pindyck (1999). 

15 For a comprehensive overview on the ADF test see Gujarati (2003), Cryer and Chan (2008), or Franke, Hardle, 
and Hafner (201 1). 
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(3.1.8) in three different specifications. First, as an autoregressive model without 
drift and deterministic time trend (AR), second, as an autoregressive model with 
drift (ARD), and, third, as a trend stationarity model with drift and deterministic 
time trend (TS). 

For each of these model specifications, the optimal lag order is determined as fol- 
lows. First, the maximum number of lags is computed according to (3.1.9) as 
proposed by Schwert (1989). Subsequently, the ADF test is programmatically run 
for all lags p = p max -, . . . , 1 and the Akaike Information Criterion (AIC) 16 of the 
regression is recorded at each step. Following Gujarati (2003), I set the lag order 
corresponding to the lowest AIC, i.e. the one that yields the best model fit. 



Pmax — 




(3.1.9) 



What follows from the consideration to allow for a higher lag order, is the question 
of the appropriate test statistic. In line with common practice, I use 



i = 



e-i 



se 



(3.1.10) 



where 9 is the Ordinary Least Squares (OLS) estimate of 0 and se its standard 
error. However, to be conservative, a lag-adjusted, unstudentised t-statistic given 
in (3.1.11) is additionally used for comparison. Here, N is the effective sample 
size, adjusted for lags. 



f 2_ *(9~1) 

§ (1-P, p p ) 



(3.1.11) 



Keeping in mind the heterogeneity of commodity price evolution in different time 
periods (see appendix A2) and the consideration of Pindyck (1999) and Gujarati 
(2003) that test results may strongly differ depending on the sample period under 
consideration, the ADF test is run on several sub-samples of the dataset and both 



16 See Akaike (1973) for further reference. 
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weekly and daily price quotes. A representative snapshot of these results is given 
in table 2. Before turning to an interpretation, let us briefly consider the choice 
of the two different sample periods analysed. First, the full sample period (1993- 
2013) is considered based on the motivation to get a preferably broad picture of 
price dynamics over a long time horizon. The choice of the second period from 
1993 to 2004 is motivated by the desire to exclude time periods that are, arguably, 
abnormal in some sense and not representative for price dynamics in general. For 
instance, given the strong link of commodity markets to global finance (Nissanke, 
2012), the financial crisis, beginning in 2007 17 , has presumably contributed to un- 
natural price dynamics thereafter. However, M. Roberts and Schlenker (2013) 
point out that a seemingly unusual ascent of commodity prices began even earlier 
in 2005, so that is appears conservative in an assessment of mean reversion to run 
the test on a period ending well ahead of these events in 2004. 

The test results suggest that a unit root can be rejected only for silver and corn and 
even in these cases, the result is not persistent across different time periods. While 
the analysis of other time periods and test statistics does not change the direction 
of these results, it must be noted that beyond the findings reported here, over the 
period from 1993 to 2006, one can reject a unit root for soybeans in a model with 
drift and on the basis of daily data for crude oil in the trend stationary model over 
the full sample period. While these findings imply little empirical support of mean- 
reverting models, Dixit and Pindyck (1994) and Pindyck (1999) would advise us 
to be cautious in the interpretation. They argue that for short time series, the ADF 
test often fails to reject a unit root in actually mean-reverting data. Additionally, it 
has to be taken into account that the failure to reject a unit root does not imply its 
acceptance. 

Hence, some further clarification is necessary before any conclusions can be drawn. 
More specifically, it is useful to analyse also the extent to which price shocks are 
temporary or permanent (Pindyck, 1999). The logic of this test follows from the 

!7 On 9 August 2007, BNP Paribas announced that it was ceasing activity in three hedge funds specialised in U.S. 
mortgage debt as they were no longer able to assess their value and risk (Elliot, 2011). 
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Table 2: 



Augmented Dickey-Fuller test results 

AR: Model without drift and deterministic trend, ARD: Model with drift only, TS: Model with drift and 
deterministic time trend. The test statistic t refers to tl, given by (3.1.10) and * ? ** ? *** indicate that a unit 
root can be rejected at the 10%, 5%, and 1% confidence level, respectively. The corresponding p-value is 
denoted by p and indicates the error probability if a unit root is rejected. Sample periods start on 02/08/1993 
and end on 30/12/2013 and 27/12/2004 depending on the sample. Data frequency is weekly. 



Resource 




AR 


ARD 


TS 


Sample 


t 


P 


t 


P 


t 


P 


Oil 


1993 -2013 


1.01 


0.92 


-1.01 


0.73 


-2.89 


0.17 




1993 - 2004 


0.48 


0.82 


-1.67 


0.44 


-2.91 


0.16 


Aluminum 


1993 -2013 


0.39 


0.79 


-2.20 


0.21 


-2.29 


0.45 




1993 - 2004 


0.75 


0.88 


-2.14 


0.23 


-2.15 


0.52 


Copper 


1993 - 2013 


1.11 


0.93 


-0.66 


0.85 


-1.60 


0.79 




1993 - 2004 


0.63 


0.85 


-0.91 


0.78 


-0.65 


0.98 


Lead 


1993 -2013 


1.18 


0.94 


-1.15 


0.68 


-1.95 


0.62 




1993 - 2004 


1.10 


0.93 


-1.30 


0.61 


-1.30 


0.89 


Nickel 


1993 -2013 


0.63 


0.85 


-1.69 


0.43 


-1.90 


0.64 




1993 - 2004 


0.97 


0.91 


-1.24 


0.63 


-1.68 


0.75 


Zinc 


1993 - 2013 


0.62 


0.85 


-1.44 


0.55 


-1.98 


0.60 




1993 - 2004 


0.38 


0.79 


-2.09 


0.26 


-2.10 


0.54 


Gold 


1993 -2013 


1.38 


0.96 


-0.01 


0.96 


-2.03 


0.58 




1993 - 2004 


0.18 


0.72 


-1.19 


0.66 


-0.81 


0.96 


Silver 


1993 -2013 


0.68 


0.86 


-0.86 


0.80 


-2.24 


0.48 




1993 - 2004 


0.06 


0.68 


-2.79* 


0.06 


-2.83 


0.19 


Wheat 


1993 -2013 


-0.67 


0.40 


-2.10 


0.25 


-2.51 


0.34 




1993 - 2004 


-0.86 


0.34 


-2.11 


0.25 


-2.16 


0.51 


Com 


1993 -2013 


-1.91* 


0.05 


-1.90 


0.34 


-2.41 


0.39 




1993 - 2004 


-1.59 


0.11 


-2.24 


0.19 


-2.64 


0.27 


Soybeans 


1993 -2013 


0.10 


0.69 


-1.33 


0.59 


-2.44 


0.37 




1993 - 2004 


-0.70 


0.39 


-2.10 


0.25 


-2.11 


0.54 


Sugar 


1993 - 2013 


0.04 


0.67 


-1.69 


0.43 


-2.09 


0.55 




1993 - 2004 


-0.39 


0.51 


-1.98 


0.30 


-2.78 


0.21 


Coffee 


1993 -2013 


0.17 


0.71 


-2.03 


0.29 


-1.97 


0.60 




1993 - 2004 


0.13 


0.70 


-1.57 


0.49 


-2.37 


0.41 


Cocoa 


1993 -2013 


0.75 


0.88 


-1.78 


0.40 


-2.54 


0.32 




1993 - 2004 


0.44 


0.81 


-2.05 


0.27 


-2.04 


0.57 
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earlier consideration of the difference in the persistence of shocks in the RW and 
AR(1) process. This time, we are not directly interested in the value of 0, but in- 
stead analyse the historical price variance to make inferences on mean reversion. 
An appropriate tool for this task are variance ratio tests of the sort proposed by 
Cochrane (1988). The test ratio can be defined as ’’the sample variance of the k- 
period differences divided by k times the sample variance of the first difference” 
(Cecchetti & Pok-sang, 1994, p. 177) and can be calculated according to (3.1.12). 



1 var{S t+k - S,) 
kvar(S t + 1 — S t ) 



(3.1.12) 



As in the ADF test, S denotes the natural logarithm of the commodity price. If a 
series is non-stationary, the numerator of this ratio should grow linearly in k and, 
hence, vr k should approach one. Conversely, if the data is stationary the k-period 
variance should approach an upper limit so that vr/. converges to zero as k grows. 
Lastly, a ratio rising above one indicates mean aversion, suggesting neither con- 
sistency with GBM nor mean-reverting processes (Pindyck, 1999). The variance 
ratio test results are reported for the previously considered two sample periods and 
for all commodities to obtain a well justified conclusion (fig. 2). 

The findings lead to a similar conclusion as the ADF test. Except for oil and silver 
none of the commodities display variance ratios that seem to tend towards zero. 
Finding consistently supporting evidence for mean reversion is even more difficult 
when considering both sampling periods together. 18 

In summary, it follows from the previous investigation that only silver exhibits 
mean reversion in both tests. Additionally, there is some evidence for mean re- 
version in corn, oil, and soybeans, however, these results are not persistent across 
tests, sample period, or data frequency. While these findings may be explained 
by a very low speed of mean reversion that can simply not be detected over two 
decades as suggested by Dixit and Pindyck (1994), it raises nonetheless questions 
about the appropriateness of mean-reverting price models in real option applica- 



18 Note that additional calculations reveal no sensitivity of the results to the choice of lag order. 
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Figure 2: Variance ratio test. This figure shows the variance ratio (y-axis) for 2...i00 lags (x-axis). 
The ratio is calculated over two time periods. First, from 02/08/1993 to 30/12/2013 (black). Second, from 
02/08/1993 to 27/12/2004. Also shown is a red reference line for a value consistent with a RW process. 
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tions. Although, it is difficult to obtain reliable estimates of typical investment 
horizons in a capital budgeting context, it seems far fetched that these are regularly 
in excess of 100 years so that mean reversion would play a role in price modelling 
according to the findings of Dixit and Pindyck (1994) and Pindyck (1999). In- 
stead, the only powerful arguments in favour of applying mean-reverting models 
are on the theoretical side. Several studies, including those of Laughton and Ja- 
coby (1993), Dixit and Pindyck (1994), Schwartz and Smith (2000), and Sick and 
Gamba (2010), propagate the use of such models mainly based on the idea that in 
an equilibrium setting, non-differentiated products, such as commodities, should 
trade at prices close to their marginal cost of production over the long-run. Indeed, 
the use of a mathematical model without theoretical support for its assumptions is 
not a convincing practice, however, I argue that it is just as little convincing to rely 
on a theory without thorough empirical support. In this regard, it is unclear at this 
stage what conclusions we should derive from the previous statistical results for 
our later task to propose meaningful stochastic commodity price models. In fact, 
before moving on, we should clarify whether the previously identified evidence 
against mean reversion can also be reconciled with economic intuition so that we 
can confidently omit the study of mean-reverting models in the remainder of this 
paper and focus on non-stationary ones instead. 

To shed some more light on the question if the assumption of stationary prices is 
sensible in a real option setting, let us first conduct a brief hands-on experiment 
to backtest the relative appropriateness of mean reversion on an intuitive level 
over the past two decades using the previously mentioned stationary Vasicek or 
Ornstein-Uhlenbeck model and GBM (see fig. 3). Since the technicalities of both 
models are not yet relevant at this stage, a detailed consideration of GBM is post- 
poned until section 4.1. and some details on the Vasicek process (as far as relevant 
for this exercise) are presented in appendix B. 

Since aluminium price data is available as of 1980, this market is a convenient 
venue for this task. To be realistic, both processes are calibrated with 15 years 
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A: Backtest 
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Figure 3: Backtesting the assumption of mean reversion. Panel A shows the historical price evolution 
of the LME aluminium price. The two shaded time periods from left to right are the model calibration 
period (02/01/1980 - 30/12/1994) and a period of high volatility (05/09/2005 to 30/08/2010), respectively. 
Also shown are volatility bounds of GBM and the Vasicek model, which represent the expected value +/- 
one standard deviation (o) . Panel B illustrates the price distribution under both processes for the date 
11/07/2008, where an intermediate price peak of US$ 3,317 occurred. The distribution is extracted from 
10,000 price simulations. 



of daily data (1980-1995) 19 and are subsequently compared to realised aluminium 
prices from 1995 until the end of 2013. Despite even higher volatility in the cali- 
bration period than thereafter, 20 it is evident that the implied distribution of prices 
under the mean-reverting Vasicek process vastly underestimates fluctuations in the 
real world. In fact, when considering the price peak on 1 1/07/2008 and the corre- 
sponding model price distributions, it stands out that observations of this kind (or 
above) have zero probability in the Vasicek model but 13 per cent under GBM. 

If we imagine now a hypothetical firm assessing the value of a natural hedge 
against high aluminium prices such as a production facility that allows to switch to 
a different, less costly input resource in times of high aluminium prices, it becomes 
clear that the inherent option to switch 21 would have been undervalued on the basis 
of assuming mean reversion in prices. In other words, if a stochastic process does 
not assign any probability to contingencies such as the price peak on 1 1/07/2008, 

19 Both models are calibrated using Maximum Likelihood Estimation (MLE). In the case of GBM this involves 
estimating the drift and volatility parameter, in the Vasicek model the long-term mean, volatility, and mean-reversion 
speed. Please see appendix B2 for details. 

20 The annualised volatility of daily log-retums was 21.6 % in the calibration period and 19.9 % in the testing period. 

21 See Tsekrekos et al. (2012) for further reference. 
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the value of flexibility, i.e. the possibility of contingent action following this event, 
is not accounted for in the real option analysis. As a result, I argue that even if 
commodity prices mean revert in the very long-run, it is the distribution of prices 
implied in mean-reverting stochastic processes that may lead to inaccurate deci- 
sion making in capital budgeting throughout time. 

Moreover, while mean reversion may have been plausible over the past 100 years 
(Pindyck, 1999), it appears to be a less coherent assumption in more recent times 
and, presumably, the future. This conjecture relies on two arguments. First, struc- 
tural supply and demand imbalances have originated from an upsurge in commod- 
ity demand from Asia and newly industrialising emerging economies as well as 
supply constraints that are expected to persist over the medium to long-term (Nis- 
sanke, 2012; M. Roberts & Schlenker, 2013). Given the resulting precarious state 
of demand and supply uncertainty, commodity prices are expected to remain highly 
dependent on shifts in market expectations and the uncertain growth prospects of 
the world economy. As a result, it is difficult to picture a realistic long-run mean 
price level in conjunction with a restricted variance of future outcomes. This argu- 
ment is supported by Tang (2012), who contends that if commodity prices are mod- 
elled as mean-reverting, the mean level needs to be stochastic. Second, commod- 
ity markets have witnessed a dramatic increase in derivatives trading and inflow 
of speculative funds since 1980, which strengthen the link of commodity prices to 
other asset classes and exacerbate the magnitude of price fluctuations (Erb & Har- 
vey, 2005; Gorton & Rouwenhorst, 2006; Nissanke, 2012). What follows is that 
commodity prices are repeatedly driven away from presumed equilibrium levels 
over extended time periods so that any meaningful assumption on mean reversion 
is difficult to hold. 

Conclusively, this section has revealed that there is little empirical support for the 
assumption of mean reversion in commodity prices over the past two decades. Fur- 
thermore, I argue that this may not be counter intuitive or at odds with economic 
intuition in the wake of structural changes in the global real economy and contin- 
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uously evolving financial markets. 



3.2 Analysis of returns 

Given the previous conclusion that the assumption of stationary prices is often dif- 
ficult to justify in today’s commodity markets, we will shift our focus to the analy- 
sis of empirical price properties relevant for non-stationary price models. The goal 
of this section is to gain an understanding if a simple non-stationary model such 
as GBM is capable of reproducing the (most relevant) statistical properties of past 
prices in future simulations. If this is not the case, it needs to be understood where 
deviations between standard model assumptions and reality occur so that solutions 
can be proposed in later sections. 

Please note that all subsequent analyses focus on the behaviour of returns instead 
of price levels. This is convenient, because the study of any time series with the 
goal to generalise the identified properties to other time periods or the parametri- 
sation of simulation models requires stationarity of data (Gujarati, 2003). While 
we have already seen that prices mostly do not fulfil this requirement, Cont (2001) 
suggests that the distributional properties of returns are more persistent over time. 
Clearly, prices and returns are just the flipside of the same quantity, however, the 
study of returns allows us to extract valuable information from the past more read- 
ily (Mantegna & Stanley, 2000). 

The remainder of this subsection is organised as follows. In the first part, some 
of the empirical properties of commodity returns are investigated. This involves 
preliminary tests on return stationarity and independence, as well as tests on the 
normality of return distributions. The second part is aimed at making a link be- 
tween the identified statistical properties of returns and stochastic price models. In 
other words, we will try to identify more explicitly how the empirical findings can 
be captured by stochastic processes. 
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3.2.1 Stylised properties 

While the list of stylised return properties one can test is long (Cont, 2001; Cont 
& Tankov, 2004), we will limit our focus to an investigation of those proper- 
ties that are directly relevant for our purpose of stochastic price modelling in the 
next section. To begin with, let us define r, as the log return at time t given by 
r, = In (P, /Pt-i). Please note that log returns are used due to their convenient mul- 
tiplicative properties and common use in financial theory (Hull, 2009). 

As a first step, we will verify whether returns are indeed stationary so that the dis- 
tributional properties of r remain similar over time. To see why this is necessary, 
imagine we want to calibrate a stochastic process to past return data and use the 
resulting parameters for future price simulations as we did in fig. 3. Clearly, for 
the results to be meaningful, the distribution of r must remain largely unchanged 
during the calibration and analysis period. Now, to test stationarity we will resort 
to the ADF test and run it on past weekly returns. As this test has already been 
discussed before, we jump directly to the results shown in table 4. In line with 
expectations, a unit root can be rejected for all commodities at the one per cent 
confidence level so that the above postulated requirement is satisfied. 

Next, it is analysed if returns are independent. This is a common assumption in 
many (non-stationary) stochastic processes and it is intimately linked to economic 
theory. In other words, if markets are efficient, we would expect current prices 
to reflect all available and predictable information. What remains to drive prices 
must be unpredictable new information that is by definition random (Fama, 1965; 
Samuelson, 1964). If this is so, we would expect the absence of any systematic 
link between today’s and tomorrow’s prices, i.e. we would expect returns to be 
independent random shocks. To verify this conjecture, we check for autocorre- 
lation in historical returns (Brigo, Dalessandro, Neugebauer, & Triki, 2007). For 
this purpose, the Ljung-Box (LB) test proposed by Ljung and Box (1978), and the 
Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) ad- 
vocated by Brigo et al. (2007) are used. 
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The LB test compares the null hypothesis that our series of returns exhibits no 
autocorrelation for a certain number of lags L with the alternative that some auto- 
correlation coefficient p (t), t = 1 , . . . ,L is different from zero. The test statistic is 
computed according to 



where T is the sample size, L the number of lags at which the autocorrelation p 
is evaluated. The test static follows a chi-square distribution with L degrees of 
freedom (Gujarati, 2003). In line with Kat and Oomen (2007), ten lags are used in 
the test. The results are shown in table 4, indicating that in the case of crude oil, 
copper, nickel, gold, and coffee, the null hypothesis of independent returns must 
be rejected. As this statistical relation between observations at different lags is 
unexpected and inconsistent with many stochastic processes, a graphical analysis 
of the ACF and PACF is provided to judge the severity of this finding. The ACF at 
lag t is given by 



where k denotes the lag at which the function is evaluated, K the maximum number 
of lags up to which the analysis is carried out, T the number of observations, and 
ju the sample mean. Before turning to an interpretation of results, let me also 
introduce the PACF as it is often considered complementary to the ACF (Brigo et 
al., 2007). While the ACF \ yields an estimate of the correlation between r(t) and 
r(t + k), PACF k indicates the correlation between these observations that is not 
explained by intermediate lags 1 1 (Gujarati, 2003). According to Brigo et 

al. (2007), the PACF is roughly a sample estimate of the quantity 




(3.2.1) 



&CF k = \ r £ (r t -ju){r t+k -ju). k = 0,1,2, ...,K (3.2.2) 

1 t— 1 




(3.2.3) 
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A: Autocorrelations of r. 



B: Partial autocorrelations of r. 
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Figure 4: Sample ACF and PACF. Panel A and B show the sample autocorrelation function and partial 
autocorrelation function for weekly copper data over the sample period from 02/08/1993 to 30/12/2013. 
Both functions are evaluated up to 20 lags in line with common practice (Brigo et al., 2007). 



where r' +1 t+k _\-J t+ k and t+k _ j denote the OLS estimates of r, and r J+J t 
given the intermediate observations r t+ 1 , . . . , i . 22 The results of this exercise 
reveal that most of the commodities, where independent returns were rejected un- 
der the LB test, exhibit one or two significant lags, suggesting that autocorrelation 
may still be acceptable (Brigo et al., 2007). An exception to this conclusion is the 
case of crude oil, where four lags enter significantly at the 95 per cent level. A 
graphical representation of the analysis is given exemplarily for the case of copper 
in fig. 4. The test results for all other commodities are available from the author, 
however, as these are not of central importance and, again, similar to those pre- 
sented for copper, they are omitted in the document. Given the results of the first 
two preliminary tests, we can draw the intermediate conclusion that the basic as- 
sumptions of stationarity and independence of returns are predominantly satisfied. 
In the next step we are interested in the distribution of returns. As a reference point, 
let us recap that the standard assumptions in GBM are log-normally distributed as- 
set prices and, accordingly, independent, normally distributed log returns (Black 
& Scholes, 1973; Merton, 1976). While this assumption has been contested in eq- 
uity markets by numerous authors, studies on commodity spot prices are rare. 23 

“The interested reader is referred to Box and Jenkins (1994) for further reference. 

23 See, for instance, Cont (2001), Schoutens (2003), Cont and Tankov (2004), Brigo et al. (2007), and Kienitz and 
Wetterau (2012) for studies on equity markets and Gorton and Rouwenhorst (2006) and Kat and Oomen (2007) for an 
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Nonetheless, all findings point in the direction that commodity returns are also 
not normally distributed. To ascertain the extent to which this notion is similarly 
supported by the dataset employed in this paper, I compute higher moments, run 
normality tests, and conduct a graphical study of empirical density in the remain- 
der of this section. 

As a first step, the sample mean (p ) and standard deviation (d) are computed for 
general reference and comparison with other studies such as those of Gorton and 
Rouwenhorst (2006) and Kat and Oomen (2007). As the calculation is standard, 
formulae are not reported here. To estimate the symmetry or skewness of returns, 
three measures are calculated. First, the most common kurtosis measure corrected 
for sample bias SI is calculated according to 

TlLih-Pf V T ( T ~ i) 



51 = 



flLiOv-#)' 



T - 2 



(3.2.4) 



where T is the total number of observations in the sample. In addition to SI, two 
complementary measures S2 and S3 are reported. According to Kat and Oomen 
(2007), these are less sensitive to outliers than SI as they measure the asymmetry 
in distribution mass based on percentiles instead of third powers. The measures 
are computed as follows: 



A 5 + As - 2P 5 o 
A 5 - As 
ft - Ao T 
lf =1 Vt - Aol 



(3.2.5) 

(3.2.6) 



Besides skewness, the fourth moment, excess kurtosis, yields an indication of the 
concentration of returns at the centre and the tails of the distribution. It is com- 
monly referred to as a measure of fat-tailedness, i.e. the extent to which extreme 
price swings occur more frequently than a normal distribution would imply (Kat & 
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Oomen, 2007). Similarly to the skewness, I report the standard measure of excess 
kurtosis adjusted for sample bias calculated as 



K 1 = 



T- 1 



(T — 2) (T — 3) 



(T + 1 ) 



fE, r =l(r,-A >) 4 






(3.2.7) 



and two less outlier sensitive metrics suggested by Kat and Oomen (2007). Again, 
these rely on percentiles of the return series to give a measure of the relative distri- 
bution mass in the tails and the centre. These measures K2 and K3 are calculated 
as follows: 



K2 = 



Psi.5 - 1*62.5 + ^* 37.5 ~ ^* 12.5 
Pl5 ~ P25 



K3 = 



P91.5 - P2.5 
Pi 5 ~ P25 



— 2.91 . 



1.23 



(3.2.8) 

(3.2.9) 



Turning to the results of these calculations reported in table 3, evidence is mixed in 
the case of skewness as some commodities exhibit positively and others negatively 
skewed distributions. Also, the signs of each of the three skewness measures are 
only consistent in six of the cases and the magnitude of skewness differs consider- 
ably. It stands out, that the highest (absolute) level of skewness is generally sug- 
gested by SI, reflecting the sensitivity of this measure to outliers (Kat & Oomen, 
2007). To examine the impact of outliers more thoroughly, I replace extreme ob- 
servations by the following replacement rule proposed by De Laurentis, Maino, 
and Molteni (2010): 



f(0 = < 



Q3 + 1.5IQR 
QI-1.5IQR 



if r(t) >Q3 + l.5IQR 
if r(t) <Ql-l.5IQR 



r(t) else. 



(3.2.10) 



The new, outlier adjusted return vector is denoted by r, Q1 and Q3 are the first and 
third quartile of r, respectively. The difference between Q3 and Q1 is the Interquar- 
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Table 3: 



Moments of commodity returns 

This table shows the first four moments of commodity returns. All indicators are calculated based on weekly 
data over the full sample period from 02/08/1993 to 30/12/2013. The sample mean and standard deviation are 
annualised with factor 52 for better comparability. All indicators are calculated based on raw data without 
any outlier treatment. 





Mean 


Std. dev 


SI 


S2 


S3 


K1 


K2 


K3 


Oil 


0.09 


0.38 


-0.55 


-0.11 


-0.10 


8.52 


0.06 


0.32 


Aluminum 


0.02 


0.21 


-0.10 


0.03 


0.01 


3.91 


0.09 


0.41 


Copper 


0.07 


0.26 


-0.62 


0.06 


-0.01 


5.57 


0.14 


0.81 


Lead 


0.09 


0.33 


-0.07 


0.07 


0.03 


5.27 


0.21 


0.91 


Nickel 


0.05 


0.36 


0.09 


0.03 


0.02 


5.32 


0.17 


0.76 


Zinc 


0.04 


0.28 


-0.46 


0.04 


-0.00 


6.15 


0.11 


0.91 


Gold 


0.05 


0.17 


0.02 


0.03 


0.01 


7.44 


0.21 


1.13 


Silver 


0.06 


0.31 


-0.97 


0.02 


-0.03 


10.16 


0.25 


1.10 


Wheat 


0.02 


0.30 


0.20 


0.06 


0.06 


5.20 


0.11 


0.70 


Com 


0.03 


0.32 


-0.48 


0.07 


0.02 


5.81 


0.17 


0.81 


Soybeans 


0.03 


0.27 


-0.68 


0.05 


-0.01 


7.54 


0.07 


0.76 


Sugar 


0.03 


0.31 


-0.20 


-0.02 


-0.02 


4.69 


0.07 


0.47 


Coffee 


0.03 


0.39 


0.60 


-0.00 


0.02 


10.47 


0.05 


0.53 


Cocoa 


0.05 


0.27 


0.16 


0.04 


0.03 


4.82 


0.24 


0.66 



tile Range (IQR). The outlier adjustment leads to an average of 39 replacements 
across the 14 commodites or, put differently, 3.7 per cent of observations per case. 
If one uses a less restrictive definition of outliers based on 3IQR instead of 1.5IQR, 
the average number of replacements per case declines to 5 or 0.5 per cent. As it 
is difficult to judge which specification is more appropriate, we will consider the 
impact of both rules where helpful. 

As suggested by Kat and Oomen (2007) also SI indicates hardly any skewness in 
the data after outlier replacement according to the 1.5IQR rule. However, using 
3IQR instead, skewness remains present particularly for measure SI. In any case, 
these findings are surprising since they contradict the popular idea of commodi- 
ties’ positive exposure to supply shocks, which would lead us to expect a generally 
positive skewness throughout all cases. 

Turning to the results for excess kurtosis, all metrics indicate values above normal. 
Even after outlier management this result is weakened, but leads to the same con- 
clusion. This is little surprising given the extreme price swings observed in past 
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data. To put the findings into context, crude oil, precious metals, coffee, and soy- 
beans exhibit values for K1 similar to the S&P 500 (K1 = 8.38). 

To investigate whether the degree of skewness and excess kurtosis just identified 
is still consistent with normally distributed returns, two normality tests are consid- 
ered next. First, the Jarque-Bera (JB) test, introduced by Jarque and Bera (1987), 
is a moment-based test exploiting the relationship between the normality of a dis- 
tribution and the two higher moments S and K. The null hypothesis that a vector of 
returns is drawn from a normal distribution is tested against the alternative that it 
is subject to a different distribution. The test statistic is given by 



and critical values are interpolated into a table generated by Monte Carlo simula- 
tion. The second test considered was introduced by Anderson and Darling (1952) 
and will be referred to as the AD test. Similarly to the JB test, it assesses the null 
hypothesis that returns are normally distributed, however, follows a different prin- 
ciple so that it serves as a useful complement to the JB test. In particular, the AD 
test derives the empirical distribution function from the return sample under in- 
vestigation and compares it to the normal cumulative distribution function. While 
there is a range of tests of this kind, 24 the AD test is found to be particularly power- 
ful and sensitive to departure from normality in the tails of the distribution (Razali, 
2011; Stephens, 1974). The test statistic is given by 



where {r\ < ■■■ < /y } are the ordered sample data points and T is the number 
of observations in the sample. The underlying principle of this test statistic uses a 
measure of total distance between the empirical and normal cumulative distribution 
functions evaluated at each of the ordered sample points. Critical values are also 

24 Alternatives include the Lilliefors, Kolmogorov-Smimov and Shapiro-Wilk test (Razali, 2011). 




(3.2.11) 
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Table 4: 



Stationarity, independence, and normality of returns 



All indicators are calculated based on weekly data over the full sample period from 02/08/1993 to 
30/12/2013. Note that outliers have not been removed from the data. 





ADF 




Ljung-Box 


Jarque-Bera 


Anderson-Darling 




t 


p 


t 


P 


t 


p 


t 


p 


Oil 


-19.00*** 


0.00 


43.20*** 


0.00 


1369*** 


0.00 


7.36*** 


0.00 


Aluminum 


-31.49*** 


0.00 


10.92 


0.36 


37.41*** 


0.00 


1.85*** 


0.00 


Copper 


-30.41*** 


0.00 


23.32*** 


0.01 


353*** 


0.00 


4.96*** 


0.00 


Lead 


-33.53*** 


0.00 


9.24 


0.51 


224*** 


0.00 


7 Q7*** 


0.00 


Nickel 


-32.29*** 


0.00 


17.49* 


0.06 


233*** 


0.00 


4.40*** 


0.00 


Zinc 


-33.68*** 


0.00 


6.49 


0.77 


454 *** 


0.00 


6.92*** 


0.00 


Gold 


-33.58*** 


0.00 


20 . 88 ** 


0.02 


851*** 


0.00 


10.16*** 


0.00 


Silver 


-33.61*** 


0.00 


15.60 


0.11 


2381*** 


0.00 


10 . 22 *** 


0.00 


Wheat 


-32.66*** 


0.00 


5.97 


0.82 


215 *** 


0.00 


4.64*** 


0.00 


Com 


-32.02*** 


0.00 


6.96 


0.73 


3gl*** 


0.00 


6.71*** 


0.00 


Soybeans 


-32.68*** 


0.00 


14.02 


0.17 


971*** 


0.00 


6.92*** 


0.00 


Sugar 


-31.92*** 


0.00 


13.19 


0.21 


129*** 


0.00 


3 03*** 


0.00 


Coffee 


-34.06*** 


0.00 


25.76*** 


0.00 


2474 *** 


0.00 


Inf*** 


0.00 


Cocoa 


-31.69*** 


0.00 


9.63 


0.47 


147*** 


0.00 


4.24*** 


0.00 



here obtained using Monte Carlo methods. Turning to the results of these tests in 
table 4, we notice that normality can be rejected at the one per cent confidence 
level under both tests and for all commodities. This result changes for the JB test 
once we account for outliers according to the replacement rule introduced above. 
Once extreme values are replaced, normality can only be rejected for crude oil, but 
none of the other commodities. However, the selection of the outlier replacement 
rule is decisive also here. If one uses the less restrictive 3IQR rule instead of 
1.5IQR, the JB test results remain unaffected compared to the case with no outlier 
management. Considering the AD test, the results remain entirely unaffected by 
outlier management regardless of the replacement rule applied. In summary, we 
are only able to accept the null hypothesis of normal returns under one of the two 
tests if a restrictive rule of so-called outlier replacement is imposed on the data, 
substituting 3.7 per cent of observations. This points us in the direction that returns 
are really not normally distributed. 

As a lest step in this question, a graphical analysis may prove helpful to gain a 
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better understanding where empirical distributions deviate from normality. In order 
to plot the historical distribution, kernel density estimation is a helpful tool that 
allows us to infer the empirical density function from independent and identically 
distributed random variables from a sample of observations (Parzen, 1962). This 
procedure is also commonly applied to asset price data 25 and, roughly speaking, 
draws a smooth line through the empirical sample frequencies observed at each 
evaluated return point (r). Accordingly, the kernel density estimate is given by 

(3 - 2 ' !3) 



where K (r) is the kernel function and h is the bandwidth or smoothing parameter. 
I follow Schoutens (2003) and use the Gaussian kernel as kernel function, which 
is given by 



K (r) = exp 




1 

V2n' 



(3.2.14) 



In line with Silverman (1986), the bandwidth is calculated using the optimal rule 
for Gaussian data given as: h = y/46 5 /3T. The results of this calculation are de- 
picted in fig. 5 for the fairly representative case of gold. In particular, the figure 
presents the empirical and normal density on a standard and log-scale. Addition- 
ally, panel C shows the empirical distribution quantiles plotted against the normal 
quantiles. In line with our previous findings, the graphical analysis verifies that in 
many cases commodity prices do not move as much as we would expect under a 
normal distribution, but at the same time we observe that the historical quantiles 
in the distribution tails are considerably larger than under the normal distribution. 
In summary, this section has shown that returns of most commodities are not nor- 
mally distributed. In some cases, positive or negative skewness is a reason for this 
observation, however, what stands out is the considerable level of excess kurtosis 
common to all commodities. While skewness may be driven by a number of ex- 



25 See Schoutens (2003) or Kienitz and Wetterau (2012). 
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A: Normal and kernel density B: Normal and kernel log density C: Quantile-quantile plot 
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Figure 5: Kernel and normal density. The plots are generated based on weekly log returns for gold spot 
prices over the sample period from 02/08/1993 to 30/12/2013. Panel A shows the kernel density (black) 
calculated according to (3.2.13) and the normal probability density (grey) with matched weekly mean (jj = 
0.0013) and standard deviation (a = 0.0359). The density is evaluated on a grid of 100 returns. Panel C 
shows a standard QQ plot, where the red line joins the first and third quartiles as a reference for normality. 



treme observations, excess kurtosis remains present even after outliers are removed 
from the dataset. As a result, the hypothesis of normally distributed returns is al- 
ways rejected under the distribution-based AD test and for less restrictive outlier 
treatment also under the JB test. 

3.2.2 Jumps and GARCH effects 

Given the findings of the previous subsection, we know that commodity returns 
are not normally distributed, most notably, due to positive excess kurtosis. In this 
section, the aim is to understand what factors cause these fat tails in return distri- 
butions so that we can come up with realistic stochastic processes that can mimic 
historical distributions. Fortunately, the list of possible ways to model fat tails and 
generally excess kurtosis in stochastic processes is not long. Accordingly, fat tails 
can be instilled into a return distribution by allowing the volatility of returns to vary 
stochastically (Bai, Russell, & Tiao, 2003) or by entering discontinuities, so-called 
jumps, into price trajectories (Brigo et ah, 2007; Schoutens, 2003). In order to see 
which of these features is likely to work behind the scenes of commodity returns, 
we will consider both possibilities in some more detail hereafter. 

Considering first the case of stochastic volatility, I compute standard deviations 
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0.08 



A: Trailing monthly volatility 



B: Daily absolute log returns 





Time Time 



Figure 6: Time-varying volatility. Panel A depicts the standard deviation of copper log returns for trail- 
ing windows of 22 days, the equivalent of approximately one trading month. Panel B, shows the absolute 
values of corresponding daily log returns. Both figures are generated with sample data from 02/08/1993 to 
30/12/2013. 



based on daily copper data for trailing windows of one month and plot the re- 
sults along with corresponding absolute weekly log returns in fig. 6. In line with 
the findings of Schoutens (2003) in equity markets, copper returns exhibit time- 
varying volatility according to panel A. Moreover, there is evidence for volatility 
clusters, i.e. a succession of periods with high or low return variance. While this 
observation is in so far interesting as it contradicts the standard assumption of time 
invariant volatility in GBM, we do not yet know to what extent it is responsible 
for the deviation from normality in the form of fat tails. To answer this question, I 
suggest to control for stochastic volatility effects and assess the resulting change in 
the normality of the distribution. In other words, if we find that returns are (more) 
normally distributed when the effect of time-varying volatility on the distribution 
is eliminated, we could conclude that stochastic volatility is indeed a relevant fac- 
tor that can capture some of the non-normality of returns through fat tails. As a 
result, it may be sensible to introduce stochastic volatility also into the price mod- 
els considered in section four. 26 

A useful tool for this kind of analysis are Generalised Autoregressive Conditional 
Heteroscedasticity (GARCH) models proposed by Engle (1982) and Bollerslev 
(1986). In particular, let us consider a GARCH(1,1) model where the commod- 

26 Note that a similar approach is used by Kat and Oomen (2007). 
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Table 5: 



Stochastic volatility and normality 

Al: Aluminium, Cc: Cocoa, Cf: Coffee, Cn: Corn, Co: Copper, Go: Gold, Le: Lead, Ni: Nickel, 01: Crude 
oil, Si: Silver, So: Soy beans, Su: Sugar, Wt: Wheat, Zi: Zinc. Results are calculated based on weekly data 
over the full sample period from 02/08/1993 to 30/12/2013. Outliers were not removed prior to calculations. 





Ol 


Al 


Co 


Le 


Ni 


Zi 


Go 


Si 


Wt 


Cn 


So 


Su 


Cf 


Cc 


DoF,_ fil 


4.75 


8.60 


5.12 


4.10 


5.03 


4.21 


3.52 


3.68 


5.11 


4.25 


4.35 


5.97 


5.07 


5.35 


DoF GARC h 


6.91 


12.57 


7.70 


7.30 


6.12 


7.21 


5.36 


5.25 


5.85 


4.87 


5.86 


7.48 


7.57 


5.49 


A AD 


0.46 


0.46 


0.43 


0.3 


0.66 


0.32 


0.41 


0.47 


0.75 


0.75 


0.53 


0.62 


0.43 


0.95 


y= a + P 


0.97 


0.99 


0.96 


1.00 


0.99 


1.00 


1.00 


0.99 


0.97 


0.92 


0.93 


0.99 


0.96 


1.00 



ity log return is represented as the following stationary variance process: r t = £,. 
where £, = o t zt and zt is a white noise process as discussed earlier in this section. 
oj is a conditional variance process given by 

aj = K + aie / 2 _i + PioLi ai + Pi<l. (3.2.15) 



This model is fitted to the historical series of returns using MLE and letting the in- 
novations follow Student’s t-distribution. Now, recap that this distribution has fat 
tails for a low number of Degrees of Freedom (DoF) and approaches the normal 
distribution if DoF is increased. If, therefore, the fitted GARCH(1,1) model im- 
plies a high number of DoF in the distribution of innovations, we can conclude that 
returns are close to normally distributed once it has been accounted for the impact 
of varying volatility. If, furthermore, DoF are lower if a Student’s t-distribution is 
fitted to the return data without controlling for stochastic volatility, we also have 
a reference for the relative importance of stochastic volatility for each commodity. 
On the contrary, if DoF do not rise at all once GARCH effects are controlled, we 
can infer that the effort to include stochastic volatility in a price model is probably 
not worthwhile. 

The results of this calculation are presented in table 5. For all commodities, the 
number of DoF in the Student’s t-distribution increases if we account for GARCH 
effects vis-a-vis the case where a stand-alone t-distribution is fitted to past returns 
without taking the effect of stochastic volatility into account. Although, this is 
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already the insight we were looking for, it is somewhat difficult to judge the rel- 
ative extent to which the non-normality of returns can be explained by GARCH 
effects so far. Recap that the t-distribution becomes asymptotically more normal 
for higher DoF in the sense that, for instance, an increase of DoF from one to 
two has a much larger impact on the similarity between the t-distribution and the 
Gaussian distribution than an increase from 99 to 100. Thus, the absolute change 
in DoF is of little use. To make the relative importance of stochastic volatility for 
each commodity more transparent, I propose A AD as a measure to quantify the 
change in total distance between the two distributions before and after GARCH 
effects are controlled. It is calculated according to (3.2.16) and shown in table 5. 

/rjiW(0,i) 0) - ^(0,1, V GA*C*) (x) I dx 
I-jF Nm (x)-F mhvfit) (x)\dx ’ C " ) 

F)v( o,i) denotes the cumulative standard normal distribution function, F t ^ Q { v garch j 
refers to the t-distribution’s cumulative distribution function with DoF retrieved 
from the fitted GARCH model, and F t ^ 0 | v /« j denotes the same distribution func- 
tion, however, the number of DoF is obtained from the best MLE fit to historical 
returns without accounting for GARCH effects. Please note that the integrals are 
numerically approximated with the trapezoid rule over the interval [-1000 1000]. 
The result of A AD will lie in the interval [0 1], where a value of zero indicates that 
stochastic volatility is entirely responsible for the non-normality of returns and one 
indicates that it has no impact on normality. If 0 < A AD < 1, stochastic volatility 
can account for some but not all the non-normality of a return distribution. 

The calculation of this measure indicates that the distance between the t-distribution 
and normal distribution declines by more then 50% in the majority of cases. Par- 
ticularly, in the case of lead and zinc stochastic volatility has a strong influence on 
return distributions, but seems to be of low importance in the case of cocoa. These 
findings are perfectly in line with those of Kat and Oomen (2007), who similarly 
find low conditional kurtosis in cocoa and persistent fat tails in a range of other 
commodities. 
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For the purpose of completeness and general interest, I also report the sum of the 
lag coefficients (y= ai + Pi) of the conditional volatility process. Similarly to our 
earlier discussion of the autoregressive coefficient in the AR(1) process, y provides 
an indication of the persistence of volatility shocks. As y has values above 0.9 
throughout, we can infer that if volatility increases (decreases), it is likely to stay 
high (low) for some time (Kat & Oomen, 2007). This is consistent with the previ- 
ous observation of volatility clusters in fig. 6. 

Having analysed the relation between stochastic volatility and excess kurtosis, it 
is natural to investigate whether there is evidence that the second source of excess 
kurtosis, jumps in prices, is additionally present in commodity return data. First, to 
gain a better understanding of this concept, I provide an illustrative price path with 
jumps in fig. 7. This path is right continuous, describing the property that if a jump 
occurs at time t, the asset price S, includes this jump. As a result, the jump pro- 
duces a discontinuity in the price trajectory of S from the left. Put differently, there 
is no chance at time t~ = t — At to unveil the uncertainty about the occurrence of 
the subsequent jump in point t. This behaviour of a price path corresponds to sud- 
den, unexpected arrivals of important information upon, which traders change their 
opinion on the value of a commodity (Merton, 1976). This type of large shock has 
to be distinguished from the many small shocks that are due to e.g. the arrival of 
less important information or minor supply and demand imbalances, which cause 
marginal changes in the asset price. In order to test for the existence of shocks, I 
apply the spectographic analysis proposed by Ait-Sahalia and Jacod (2012). While 
this test is devised for the application to high frequency data, I show graphically 
that it is also useful to gain some basic insights about the empirical jumpiness of 
commodity prices measured at a lower frequency. Although, a detailed explana- 
tion of the test method is beyond the scope of this paper, let us briefly consider the 
general principle. 

For any historical series of n asset prices (5) over the interval [0 T], we can 
observe discrete price changes AS'(f) = S(t) — S(t — 1) for each observation t. 
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Figure 7: Discontinuous price trajectory. This graph depicts a price path generated with the jump diffusion 
model proposed by Merton (1976) based on the following parameters: Initial asset price = 100, risk free 
interest rate = 0.05, time horizon = 1 year, volatility = 0.05, volatility of jumps = 0.05, mean of jumps = 0, 
jump intensity =10, discretisation steps = 250. A more detailed explanation of this model is given in section 
4.1.3 of this paper. 



These are to be contrasted with the actual, unobservable jumps of S denoted as 
AS = S(t) — S(t)~, where S(t)~ is the asset price one infinitesimally small quantity 
of time before the jump occurs. Now, we assume that the stochastic data generat- 
ing process of p is an Ito semimartingale. 27 A semimartingale can be decomposed 
into the sum of a drift, a continuous Brownian-motion part, and a discontinuous 
jump part. More formally we can write 



S t — .S'o T- f b s ds+ f Oyf/VVj - 

Jo Jo 

drift Brownian part 

+ / / xp(ds,dx ), 

Jo J{|jt|>e} 

S v / 

big jumps 



77 

Jo J{\ 



{W<e} 



x(p — v)(ds,dx) (3.2.17) 



small jumps 



where p is the jump measure of S and its predictable compensator is the Levy 
measure v. 28 The distinction between small and big jumps ( x ) according to the 
threshold e is arbitrary, but fixed. To learn which of the different components of 
the semimartingale are necessary to produce the historical returns and what their 

27 See Cont and Tankov (2004) or Bjork (2009) for further reference. 

28 The jump compensator ensures that the process remains a martingale in the presence of jumps (Battauz & Ortu, 
2009; Bjork, 2009). A more detailed discussion on martingales and risk- neutrality is given in section four. 
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relative importance is, the following measurement device is used: 

T/A„ 

B(k, Un , An) = £ |A?S|*1 { | A?S |< m „ } . (3.2.18) 

1=1 

The quantity B is computed for different powers ( k ), truncation levels (m), and sam- 
pling frequencies (A). For instance, the authors explain that powers of k < 2, will 
put more emphasis on the continuous component of the process, while k > 2 will 
accentuate the jump parts. Truncating the increments at u n to exclude e.g. large or 
small increments ((3.2.18) is changed accordingly) enables us to extract informa- 
tion on different types of jumps. Besides, adjustments in sampled return frequen- 
cies A„ help to distinguish whether variations converge to a finite limit, to zero, or 
diverge to infinity so that inferences can be made about which component of the 
semimartingale dominates at a particular power. Based on these different compu- 
tations, a number of test statistics can be computed, which indicate the presence of 
jumps (Sj), the degree of jump activity (Sff), the presence of continuous Brownian 
motion ( S w ), and the relative magnitude of variation in S that is attributable to a 
jump or Brownian component ( QVs p iu ). 29 As the test is originally developed for 
high frequency data, whereas we are interested in price behaviour on a daily (or 
lower) frequency, I decide to show the test statistics not only for daily historical 
commodity returns but also for a corresponding jump diffusion price path similar 
to the one presented in fig. 7. As we know the properties of the jump diffusion 
path, it is clear that the test should expose both, the jump and diffusion compo- 
nent, so that we get a better feeling for the overall reliability of the results in our 
low frequency context. Fig. 8 shows the results for both, the jump diffusion test 
path and the exemplary and representative case of coffee. Please note that, again, 
the other test results are available from the author, however, since little additional 
insights would result and in order not to exaggerate the length of this document, a 
presentation is omitted. Turning to the findings, we may notice first that the test 

29 Please see Ait-Sahalia and Jacod (2012) for further explanation. Implementation details and sample code are 
available from http: / /www. princeton .edu/~yacine/research. 
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Figure 8: Spectographic jump test. The four columns of the plot refer to the four test statistics of the 
spectographic jump test. Sj: Presence of jumps, Sfa- Degree of jump activity, S w : Presence of a Brownian 
motion, QVsplit '■ Relative importance of jump and Brownian component. The top four plots refer to the 
Merton jump diffusion model as a benchmark and reference point. The bottom four plots refer to the daily 
coffee spot price sampled over the period: 02/08/1993 - 30/12/2013. Note that daily data is used instead 
of weekly prices as the interpretation of jumps at such low frequency is problematic (Ait-Sahalia & Jacod, 
2012 ). 

results are in line with expectations for the test path. In particular, the presence 
of jumps is confirmed, there is evidence for infinite activity jumps (small Brow- 
nian increments), finite activity jumps (irregular large jumps from the compound 
Poisson process) 30 , the presence of a Brownian motion or diffusion component is 
confirmed, and, lastly, we see that the variations in S are caused by both a jump 
and Brownian process. In the case of coffee, we may notice some evidence for the 
presence of jumps, however, the relative importance of the Brownian component is 
clearly higher than under the stylised model path. As these findings are almost per- 
fectly representative for all other commodities in the sample, we can conclude that 
jumps, although not as clearly as stochastic volatility, may need to be considered in 
commodity price modelling to mimic the fat tails of historical return distributions. 



30 As mentioned previously, a more detailed explanation of the Merton model is postponed to section four. For 
further reference see Merton (1976). 
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4 Modelling commodity prices 

On the basis of the previous statistical findings, it is the objective of this section 
to propose models for the stochastic simulation of commodity prices that are the- 
oretically consistent with empirical data. In this regard, recap that mean reversion 
cannot generally be proven in past prices over the last two decades. Even if a much 
longer time horizon would allow us to reject a unit root in prices, we concluded 
that the statistical properties of mean reversion seem to neither fit reality over the 
recent past (as shown by the backtest in fig. 3) nor the future outlook of uncer- 
tain and persistently volatile commodity markets. It has also been identified that 
past returns are not normally distributed, most notably, due to excess kurtosis and 
that time-varying volatility and jumps in price trajectories are both possible causes 
for non-normality. As discussed in section one of this thesis, the most common 
stochastic processes in real options applications assume either mean reversion of 
prices or normally distributed returns. Hence, we will assess the merits of alter- 
native price processes that have already found widespread application in pricing 
equity and fixed income derivatives also in the case of commodity prices. In this 
sense, I follow Brooks and Prokopczuk (2013), who similarly draw on the often 
more developed and broad literature on equity derivatives pricing and transfer con- 
cepts such as stochastic volatility models to commodity markets. The remainder 
of this section is organised as follows. In the first half, the popular case of GBM 
is considered as a reference point and convenient example to clarify some funda- 
mental concepts of stochastic calculus and risk-neutral simulation in Monte Carlo- 
based valuation. Thereafter, we will step-by-step relax some of the restrictive as- 
sumptions of GBM on the normality of returns, constant volatility, and continuous 
price trajectories by introducing theoretically more realistic stochastic processes 
in the following subsections. In the second half of this section, each process is 
evaluated in terms of goodness of fit to past commodity time series so that first 
conclusions can be drawn on the relative appropriateness of different models. 
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4. 1 Stochastic processes 



The idea to describe asset prices in a speculative market using a stochastic pro- 
cess 31 was first proposed by Bachelier (1900), but gained little attention until the 
1950s when empirical studies revealed that serial correlation of asset price returns 
was negligible (Mantegna & Stanley, 2000) and the Efficient Market Hypothe- 
sis (EMH) was postulated by Samuelson (1964). Since then, various stochastic 
processes have been suggested for different applications. Simple examples of a 
stochastic process include, for instance, the RW and AR process as well as GBM, 
the Vasicek, or the Merton model we have already touched upon peripherally in 
the context of empirical data analysis. In this part of the paper, we will focus more 
explicitly on stochastic processes for commodity price simulation. To begin with, 
a more detailed analysis of GBM is given - a process that the reader will be fa- 
miliar with already. The reason for this choice of a starting point is twofold. First, 
it has already been mentioned that GBM is a popular model in the context of real 
options applications so that the relative performance of this process to others may 
be of general interest. Second, GBM lends itself to the preliminary study of some 
central concepts in stochastic calculus and Monte Carlo simulation without exag- 
gerated complexity. In particular, it is useful at this stage to elaborate on three top- 
ics in more detail upfront. First, the continuous time representation of a stochastic 
process in the form of its Stochastic Differential Equation (SDE). Second, how to 
solve the SDE via stochastic calculus in order to simulate the process on a discrete 
time grid of days, weeks, or months. Third, the difference between process sim- 
ulation under the historical probability measure P and the risk-neutral probability 
measure Q, which is later required for the valuation of contingent claims. 

Let us first consider the continuous time representation of GBM. Accordingly, the 



31 A stochastic process can be defined as a sequence of chronologically ordered random variables {Put > 0} (Franke 
et al., 201 1). See also Schoutens (2003) and Bjork (2009) for additional reference. 
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asset price S is governed by the expression 



S(t)=S( 0)+ fs(s)pds + f S(s)odW p (s), (4.1.1) 

Jo Jo 

where /j refers to the drift or instantaneous rate of return of S, a denotes (constant) 
volatility, and dW p is a standard Brownian motion under the historical or real 
world probability measure P with respect to the information structure or filtration 
F. In other words, dW p can be interpreted as the previously discussed white noise 
process or simply a normally distributed disturbance term that is independent of 
everything that has happened up to time f. 32 While this continuous time repre- 
sentation is an elegant building block for derivative pricing formulae such as the 
Black-Scholes model, we rarely encounter real options that admit closed form, 
continuous time solutions. As a result, we either need to approximate GBM in 
a lattice framework or rely on Monte Carlo simulation of the adequately discre- 
tised process. 33 As the Monte Carlo-based approach will prove most flexible for 
our purposes of comparing different processes, for which no lattice representation 
may be available, we need to discretise (4.1.1). This brings us to the second main 
question we want to answer in this introductory section, namely, how to solve the 
SDE for a discrete price path. 

Many stochastic processes do not admit an exact closed form solution to their SDE 
and need to be solved via some approximative discretisation scheme. 34 This is not 
the case with GBM. However, the presence of an integral with respect to dW p (s) 
prevents us, from using the usual rules of differential calculus (Battauz & Ortu, 
2009). Instead, we can first write (4.1.1) in short notation as: 

dS(t) =/jS(t)dt + aS(t)dW{t). (4.1.2) 



32 See Bjork (2009) for the formal requirements that must hold for W to be a Brownian motion. 

33 See Cox et al. (1979) for technical reference on lattice approximation, Trigeorgis (1996) or Copeland and An- 
tikarov (2003) for lattice-based real option valuation, and Tsekrekos et al. (2012) or Yungchih Wang (2011) for the 
application of Monte Carlo simulation in investment valuation. 

34 Common schemes include the Euler, Predictor-Corrector, or Milstein scheme. See Kienitz and Wetterau (2012) 
for further reference. 
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According to Bjork (2009), we can apply Ito’s formula to this type of stochastic 
differential. Accordingly we need to compute 



d f(t, s (t)) - | ^ + 1 -a 2 ^dt + o^dW F . (4.1.3) 

In line with the assumptions of GBM we set / (f , 5) = In 5 and compute the respec- 
tive differentials as: 

5 / 5f 1 5 2 / = _ 1 

5 1 ’ 55 5’ 5.S' 2 .S' 2 ' 



Substituting these into (4.1.3) and reversely substituting /(f,5) = In 5 yields: 



In 



S(0 

5(f-l) 



1 1 



l -^2 ) 5 ' 



1 



(4.1.4) 



5 2 

Exponentiating and writing the result in terms of t as opposed to dt, we get 
5(f) = 5(0) exp ^o 2 j f + oW p (f) 



(4.1.5) 



which is the familiar GBM representation in discrete time. Although, this formula- 
tion allows to simulate stochastic price trajectories under the real world probability 
measure, it does not yet allow us to price derivatives in a Monte Carlo framework. 
This leads to the third question of how to simulate prices in the risk-neutral Monte 
Carlo-based valuation of contingent claims. 

Intuitively, we know from the findings of Cox et al. (1979) that the possibility of 
creating a replicating portfolio allows the option writer to eliminate all risk arising 
from the uncertain movement of the underlying asset, on which the derivative is 
written. While this logic allows to pin down the value of the option, it also implies 
that in a risk-neutral world all assets would earn the risk-free rate of interest since 
risk-neutral investors would not demand a higher rate of return for holding risky 
assets (Glasserman, 2003). As a result, it is intuitive that contingent claims valua- 
tion using the Monte Carlo method will require us to simulate asset prices in such 
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a way that the simulated asset return corresponds in expectation to the risk-free 
interest rate. More formally, we require an equivalent probability measure Q such 
that under this new measure all discounted risky asset prices become martingales 
in expectation ’ 5 such that 



E” 



S(/ + l)l S(t) 

B(r + l)J “ B(ty 



(4.1.6) 



where B(t) = exp (rf) is a continuously compounded bond with a value of one 
and r the risk-free rate of interest (Battauz & Ortu, 2009). Roughly speaking, we 
are trying to adjust the measure, under which we formulate expectations on future 
asset prices in such a way that an otherwise risky security yields a return equal 
to the risk-free rate of interest in expectation. By the Girsanov Theorem, it can 
be shown that the existence of such a probability measure Q is equivalent to the 
existence of a market price of risk process given by: v = (p — r) /o. 36 Now, with 
reference to (4.1.5), it should hold that IV® = W + vt. To obtain the risk-neutral 
representation of the process, we can substitute fj = r + <jv and collect terms in the 
exponent of (4.1.5) to obtain 

^r + ov — t + adW p (t)= ^r— ^o 2 ^ t + o (dW p (t) + Vt^j 

= - ^cr ^ t + odW Q (t) . 



Replacing the exponent of (4.1.5) with this result, we see that our earlier intuition 
that an asset should yield the risk-free rate in the risk-neutral setting is confirmed. 37 
To put this result into context, we can value a typical European call option deriva- 
tive by simulating an appropriate number of price paths (or only the prices rele- 
vant to compute the payoffs) under the risk-neutral measure, calculate the terminal 



35 See, for instance, Mantegna and Stanley (2000) or Bjork (2009) for further reference on martingale theory. 
36 Note that under the existence of this process, the usual no-arbitrage assumption holds and the market is said to be 
complete if this process can be uniquely determined (Schoutens, 2003). 

37 Please note that in the remainder of this section, the superscripts for the historical and risk-neutral probability 
measure are omitted for notational convenience. 
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payoffs, discount them to the present at the risk-free interest rate, and obtain a val- 
uation by taking the average across paths (Battauz & Ortu, 2009). Note that in the 
presence of a convenience yield or storage costs, these have to be subtracted from 
the risk-free interest rate in the drift term (Glasserman, 20 03). 38 What remains to 
be considered, according to the logic adopted in the next subsections, is a sensitiv- 
ity analysis of asset price paths and the return distribution with respect to changes 
in the input parameters /j and a. Since, however, this exercise would add little new 
insights in the case of GBM, the analysis is omitted. 

In summary, note that while it took us some time to go through the solution of the 
SDE and the concept of risk-neutral simulation for a Monte Carlo setting, these 
considerations are to a large extent valid for the more complex stochastic pro- 
cesses presented next. Even if exact solutions to the SDEs are not available in 
some cases and require the use of some discretisation scheme, the idea of solving 
the SDE to work in discrete time remains unchanged. This logic applies equally to 
the martingale correction, which is, in many cases, more difficult to derive than for 
GBM, but is intuitively equivalent. As a result, I will omit lengthy derivations in 
the following three subsections to focus more on intuition, application, and those 
formulae directly relevant to obtain final results. 

4.1.1 Stochastic volatility 

Earlier, we concluded that the volatility of commodity returns is not constant and 
that it accounts for fat tails in the distribution of returns. As Pindyck (1999) sug- 
gests that ’’the GBM assumption will be appropriate only if [...] volatility is rela- 
tively constant” (p. 24), it seems worthwhile to allow volatility to follow a stochas- 
tic process itself. A popular model for this purpose is the one proposed by Heston 
(1993), which has already been successfully applied to commodity markets by 
Brooks and Prokopczuk (2013). The model is characterised by the following set 



38 For the possibility of return shortfall, applying in some cases to non-traded capital investments, please refer to 
Tsekrekos et al. (2012) for further reference. 
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of SDEs: 



dS(t) = vS ( t)dt + VW)S (t) dWi (t) 
dV (t)=K(&-V(t))dt+v^/vJtjdW 2 (t) 



(4.1.7) 



{dW u dW 2 ) = p dt. 

In line with previous notation, S is the asset price, V the variance process, 0 the 
long-term variance, v the volatility of variance, and p the correlation between the 
Brownian motions W\ and W 2 , which allow us to model the commonly observed 
relationship between volatility and asset prices. 39 

In order to simulate the Heston process in discrete time, we cannot solve the SDE 
by stochastic calculus as we did in the case of GBM. Instead, a number of sam- 
pling schemes have been proposed. As the exact scheme introduced by Broadie 
and Kaya (2006) is computationally expensive, I use the Quadratic Exponential 
scheme of Andersen (2008) as described by Kienitz and Wetterau (2012). Since 
the analytic representation of the QE scheme is somewhat lengthy and of little in- 
tuitive interest, please refer to appendix B4 for details. 

More interesting in our context, however, is the probability density of returns im- 
plied by the model. This will be used for later model calibration and helps us to 
understand exactly how the value of each input parameter affects the shape of the 
return distribution. As the probability density is not given explicitly, we have to 
compute it from the model’s characteristic function via the fast Fourier transform 
algorithm (see appendix B3). The characteristic function is obtained from Kienitz 
and Wetterau (2012) and can be written as follows: 40 

<p Heston (u,t) = exp (A H ( u,t ) +B C ( u,t)V (t) + iuX (t)) (4.1.8) 



39 See, for instance, Schoutens (2003) for the case of equity markets and Richter and Sorensen (2003) for commodi- 
ties. 

40 For further reference on characteristic functions in the context of stochastic processes see Schoutens (2003) and 
Cont and Tankov (2004). 
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A H (u,t) = 
B a (u,t) = 
G = 
D = 



K0 



(k — pvui — D) t — 2 log 



Gexp (— Dt ) 
G- 1 



- 1 



k - pvui - D 
CO 2 

K - pvui - D 



1 — exp(— Dt) 



1 — Gexp (—Dt) 



K-p vui + D 

\J (K-p vui) 2 + u (i + u)v 2 . 



With respect to notation, u refers to an appropriate grid of values, for which the 
function is evaluated and i is the imaginary unit from complex numbers theory. 
The relation between model inputs, the distribution of returns, and sample paths for 
the price process are shown in fig. 9. As we would expect, higher levels of 0 lead to 
a wider distribution of returns, a more volatile asset price path, and, by definition, 
a higher level of long-term variance. Considering the mean reversion speed of the 
variance to its long-term level (k), we see first that we obtain a peaked distribution 
for low values. This is intuitive since we have set 0 greater than the initial variance. 
For a low K, the variance remains low for a longer time, leading to more returns 
clustered around the centre of the distribution. This logic is confirmed by the asset 
price and variance paths. The interpretation of the volatility of variance (v) is less 
straightforward, but offers, in any case, another lever to control the overall level 
of volatility (Kienitz & Wetterau, 2012). Eventually, the correlation between asset 
prices and volatility (p) allows to control the skewness of the return distribution. 



4.1.2 Jump diffusion 

Following from the study of jumps in commodity prices, another potentially fruit- 
ful direction to pursue is the use of diffusion processes similar to GBM, but aug- 
mented by a jump component. In this section, we will consider two processes of 
this kind. First, the classical jump diffusion process by Merton (1976), which we 
have already used as a benchmark to judge the reliability of the test for the presence 
of jumps in past prices. Second, the model by Bates (1996), which can be seen as 
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Figure 9: The Heston model: Parameter sensitivities. T = 2, S(0) = 100, V(0) = 0.01, 0 = 
0.005,0.05,0.1, K = 0.001,0.1,0.5, v = 0.01,0.1,0.25, p = —0.9, 0,0.9. The base case refers to the cen- 
tral value in the previously specified ranges. If the sensitivity of output is evaluated with respect to one of 
the (four) variables under consideration, all other parameters will be set to the base case value. Note that all 
paths are calculated based on the same pseudo random number stream and are calculated under the physical 
probability measure P. Densities are evaluated on a grid of 1000 points. 



a combination of the Heston and Merton model allowing for stochastic volatility 
and jumps. 

Under the Merton model, price dynamics are governed by the following set of 
equations: 



dS (f ) = /jS (t) dt + oS (t) dW ( t ) + S (t) {Y - 1) dN ( t ) (4.1.9) 

F= W exp^-^ + OyzV z~5\£(0,1). 

The parameters fj and a correspond to the drift and volatility as in the Black- 
Scholes model. Y governs the log-normal jump size distribution, with mean jump 
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size fjj and volatility of jump size oj. N(t) is a Poisson process, determining 
the random occurrences of jumps with intensity X. By Ito’s formula for jump 
diffusions, we can solve for the discrete representation of the process in a similar 
way as for GBM and get 



yields a martingale to simulate the process under risk-neutral dynamics for the 
application in Monte Carlo-based valuation. 

As in the Heston model, we have to compute the return density for the Merton 
process from its characteristic function, which is given by: 



As we would expect, this is just the simplified product of the characteristic function 
from the normal distribution (for the Black-Scholes or GBM part of the process) 
and the logarithmic jump part. To see now the impact of different parameters on 
the return distribution and price paths, let us turn to fig. 9. As in the case of GBM, 
we skip the consideration of parameters /j and o. The jump mean controls the 
skewness of the distribution. This is sensible as we would, for instance, expect the 
return distribution to show more probability mass around the left tail if the jumps 
have on average negative values. The volatility of jumps controls the peakedness 
of the distribution. Little volatility results in few extreme observations and a con- 
centration of probability around the centre of the distribution. Finally, the effect of 
jump intensity is similar to jump volatility. If fewer jumps occur, the number of 
extreme observations decreases so that returns cluster around the centre. 




(4.1.10) 



Kienitz and Wetterau (2012) show that setting the drift to 




(4.1.11) 




(4.1.12) 
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Q. 






Figure 10: The Merton model: Parameter sensitivities. T = 10, 5(0) = 100, /j = 0.1, g = 0.1, /jj = 
—0.15,0,0.15, Cj = 0.001,0.05,0.1, X = 0.01,0.5,2. See footnote of fig. 9 for more details. 



We have seen that the incorporation of volatility or jumps in a stochastic process 
offers more control over the shape of return distributions and the dynamics of price 
paths. Thus, it may be useful to take this idea one step further. In particular, the 
Bates model can be seen as a Heston model with the addition of log-normal jumps 
from the Merton process. While the process was originally motivated by the de- 
sire to model smiles in the implied volatility surface even for short-dated financial 
options (Kienitz & Wetterau, 2012), we will investigate whether the additional pa- 
rameters also lead to a better model fit to past return data. According to Bates 
(1996), the SDE governing the evolution of this process is given by 

dS(t) =/jS{t)dt + ^/VJt)S(t)dWi{t) + {Y-l)S(t)dN(t) (4.1.13) 
dV(t)=K{&-V(t))dt+v^/vJt)dW 2 (t) 

(dW \ . dW 2 )=pdt. 

With respect to notation, we notice immediately the similarity to the Heston pro- 
cess and the additional jump part inherited from the Merton model. For the dis- 
cretisation of the Bates model, I use the Quadratic Exponential scheme known 
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from the Heston model but, augmented for jumps as described by Kienitz and 
Wetterau (2012). What remains to be discussed is the distribution of returns. To 
this end, Kienitz and Wetterau (2012) show that the characteristic function of the 
Bates process is given by 



q>Bates = cp Heston • exp (kt (-pjiu + ex p (A) - 1)) (4.1.14) 

1 , 

A = iulog{ 1 +pj) + — Gjiu ( iu — 1 ) . 

The first part of (4.1.14) is just the characteristic function from the Heston model 
and the second part the characteristic function of the log-normal jump part from 
the Merton model we have seen earlier in simplified form. Since we have already 
analysed the sensitivities of all parameters of the Bates model separately for the 
Heston and Merton process, this part is skipped here to avoid duplication of work. 

4.1.3 Levy processes 

One common characteristic of the previously considered processes is the presence 
of a Brownian motion W which affects the price process continuously. Even when 
jumps were added, we only introduced some irregular discontinuities in otherwise 
continuous price trajectories. Another class of models that does not necessarily 
contain a Brownian part is known as exponential Levy models, in honour of Paul 
Levy, the pioneer of the concept (Schoutens, 2003). These processes are pure jump 
models based on infinitely divisible Levy distributions instead of the Gaussian 
normal distribution. Similarly to Brownian motion based processes. Levy models 
have independent and stationary increments, but are more flexible than models 
based on the normal distribution in the way they can account for skewness and 
kurtosis of returns (Cont & Tankov, 2004). For this reason, it may be worthwhile 
to study the applicability of such models in commodity markets and real option 
valuation. While the array of popular Levy processes is long, 41 we will study the 



41 See Schoutens (2003) and the sources cited therein for a more comprehensive overview of Levy processes. 
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Variance Gamma (VG) process by Madan and Seneta (1990) and Boyarchenko and 
Levendorskii (as cited in Kienitz & Wetterau, 2012), the Normal Inverse Gaussian 
(NIG) process, and the NIG process with stochastic clock governed by a Cox- 
Ingersoll-Ross (CIR) process as outlined by Schoutens (2003) and Kienitz and 
Wetterau (2012). 42 

Beginning with the VG process, we can choose between two representations. First, 
the process can be expressed as the difference between two Gamma processes U ( t ) 
and D ( t ) so that S ( t ) := U (t) — D ( t ). Alternatively, we can follow Brigo et al. 
(2007) and represent the VG process as a time-changed Brownian motion. 43 For 
notational convenience, we pursue the latter of the two possibilities and write 



where the difference to the earlier considered SDE of GBM lies in the term g (f), 
which characterises market activity time (Brigo et al., 2007). More specifically, 
market time is assumed to follow an increasing random process with stationary 
increments that reconciles with calendar time on average over a given time span. 
This process is driven by random variables drawn from a gamma distribution as 
follows: {g(f)} ~ r(f/v,v). Given the SDE of the VG process and knowing the 
solution to the GBM SDE, it is intuitive that the solution to (4.1.15) in discrete 
time and under the historical probability measure is given by 



In order to obtain risk-neutrality, Kienitz and Wetterau (2012) show that we need 
to replace the historical drift measure /j by the expression 



S ( t ) = tiS (t)dt + Qdg(t) +odW ( g(t )) , 



(4.1.15) 



S(t) =S( 0)exp ^idt + Qg(t) +Oy/g(t)W (f)] . (4.1.16) 




(4.1.17) 



42 See Cox, Ingersoll, and Ross (1985) for details on the stand-alone CIR process. 
43 For further reference on Brownian subordination see Cont and Tankov (2004). 
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where r denotes the risk-free interest rate as usual. For the probability density of 
returns, the VG process admits an explicit solution, taken from Brigo et al. (2007). 
Accordingly we can calculate 



fvG (■*) 



2e 



9 (x-y) 



o v / 2rcv^r(i) 



1*-H 

V / ^ + o 2 



K, 




(4.1.18) 



where x refers to a grid of returns, where the density is evaluated, and (•) is a 
modified Bessel function of the third kind with index r| specified in the subscript 
of K in (4.1.18). To understand how the input parameters affect returns and price 
trajectories, fig. 1 1 provides the usual overview. Note that according to Brigo et 
al. (2007), /j and o play an identical role in the VG process and GBM so that 
we can omit the analysis of these parameters also here. The results show that 







Figure 11: The VG model: Parameter sensitivities. T = 10, 5(0) = 100, n = 0, o = 0.15, 0 = —0.1, 0,0.1, 
v = 0.05,0.3, 1. See footnote of fig. 9 for more details. 



the parameter 0 gives us control over the skewness and v influences the peak of 
the distribution. As expected, we observe rising price paths for positively skewed 
returns (here, more positive returns occur) and less volatile price path for a highly 
peaked distribution (fewer extreme events occur). 
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Another popular and flexible Levy process is the NIG model. Similarly to the VG 
process, the model can be expressed as a time-changed Brownian motion, where 
the time process is governed by an Inverse Gaussian distribution. Letting S(t) 
denote the price path, then, for time s < t and o.ieR, we can write S (t) -SM~ 
ig (a(t — s,b)) (Kienitz & Wetterau, 2012). As the SDE representation follows 
the same logic as under the VG process, we can, therefore, confidently move on to 
the discrete time representation without any loss of additional insights. The asset 
price under the historical probability measure at time t is governed by the following 
expression: 

S(t) =S (0)exp \pt + P5 2 T + SVYW (r)j (4.1.19) 

T = 0 + J 2 ( e - 2 - ^29 t 2 Z + 0V) , z~5\C(0,l) 

e = — 1 

hy^-p 2 



where a, p, and 5 are the NIG specific input parameters. Based on Kienitz and 
Wetterau (2012) the martingale adjustment is given by 



;U = r + 5^a 2 -(p + l) 2 -^a 2 -p 2 ^ . 
and the return density can be computed explicitly according to 



(4.1.20) 



f NIG M 




(4.1.21) 

where K is, again, the modified Bessel function of the third kind with an index of 
one. 

Next, an overview is given about how the three NIG specific parameters affect 
the distribution of returns and price paths. The first parameter a determines the 
variance of the distribution. Low values lead to high volatility and fat tails, which 
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Figure 12: The NIG model: Parameter sensitivities. T = 10, 5(0) = 100, /j = 0, a = 2,5, 10, p = 
—4, 0.5, 4, 8 = 0.05,0.2,0.8. See footnote of fig. 9 for more details. 



translate into extreme price movements. (3 affects all moments of the distribution. 
In our example, it can be observed that high values for (3 raise the expectation and 
skewness while the opposite is true for low betas. The third parameter 8 controls 
the peakedness and variance of the density, such that low parameter values lead to 
more probability mass in the centre and high values in the tails. The influence on 
the volatility of price trajectories is clearly observable. 

The next and last process that is considered is in analogy to the Bates model, where 
we were able to summarise the characteristics of jumps and stochastic volatility in 
one process. Here, we want to augment the NIG (jump) process with stochas- 
tic volatility. This is not only consistent with our earlier observation of time- 
dependent volatility in commodity returns, but also fits the intuitive argument that 
the environment changes over time so that a static level of uncertainty is difficult 
to justify (Schoutens, 2003). To introduce stochastic volatility into the NIG pro- 
cess, we employ the logic developed by Carr, Geman, Madan, and Yor (2003), 
who contend that ’’random changes in volatility can [...] be captured by random 
changes in time” (p. 351). In other words, the level of uncertainty can be increased 
(decreased) by speeding up (slowing down) the rate at which time passes. Without 
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low base high 

Figure 13: The NIG-CIR model: Parameter sensitivities. T = 10, 5(0) = 100, /j = 0, a = 5, p = 0.5, 
8 = 0.2, X = 0.01,0.5,3, k = 0.005, 1,3, r\ = 0.05,0.5, 1. See footnote of fig. 9 for more details. 



digressing to much into technicalities, we want to compute 

S{t)=S (0) exp (/ ut +L(Y (r))) , (4.1.22) 

where L is the NIG process and Y is a monotonically increasing and mean-reverting 
square root CIR process (Carr et al., 2003; Kienitz & Wetterau, 2012), which gov- 
erns the stochastic evolution of time. Note that it must be monotonically increas- 
ing as time cannot go backward. As the analytical solutions to the formula stated 
above and the return density are somewhat cumbersome, these are not presented 
here for more convenient reading. However, I provide details on the implementa- 
tion in Matlab file pathjiigcir and kindly refer the reader to Carr et al. (2003) and 
Kienitz and Wetterau (2012) for further reference. Instead, let us directly proceed 
with the graphical analysis of parameter sensitivities. Please note that since we are 
already familiar with the parameters a, (3, and 8 from the standalone NIG process, 
we limit our attention only to the parameters governing the stochastic time change 
or volatility. 

As we can observe, all of the three parameters influence the kurtosis of the distribu- 
tion, however, in a different way and with diverse impacts on the price trajectories. 
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While it is generally in line with expectations that stochastic volatility affects the 
”fat-taildness” of the distribution, a distinctive interpretation of each parameter is 
less intuitive. For our purposes, however, it suffices to notice that we gain even 
more control over the distribution shape at the cost of a additional parameters with 
similar influence on the distribution so that calibration results could be unstable. 

4.2 Model selection 

Now that have a range of models at hand, the next step is to assess their relative 
goodness of fit to past commodity return data. To do so, this section is split into two 
parts. First, a flexible and effective calibration method is introduced to obtain the 
optimal parameter sets for each of the previous stochastic processes across the 14 
commodities in the sample. Second, given these parameters we check the relative 
ability of each process to capture the statistical properties of past commodity return 
data based on an in-sample and out-of-sample test to detect possible overfitting. 

4.2.1 Calibration 

The calibration of advanced stochastic processes in our context is in some sense 
not a trivial task. As the previously discussed models are mostly envisaged for the 
pricing of market traded options, existing literature propagates the calibration of 
these processes to a basket of such options with different strike prices and maturi- 
ties to capture information about the implied volatility surface. In particular, it is 
common to algorithmically minimise some measure of distance (e.g. the root mean 
squared error) between option prices in the basket and model prices subject to the 
input parameters. This procedure directly yields the risk-neutral parametrisation 
of the process, based upon which the pricing model is build. While the calibration 
to option prices has the advantage to be forward looking in the sense that prices 
reflect the aggregated opinion of all market participants on future payoffs, data 
availability on a basket of commodity options may resemble an obstacle for cor- 
porate finance practitioners in non-financial institutions. As a result, I presume 
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that in a real options context it is much more convenient and flexible to be able 
to back out the parameters from a series of historical commodity prices directly. 
Along these lines a number of methods have been suggested. Gibson and Schwartz 
(1990) propose a seemingly unrelated regression model to fit mean-reverting mod- 
els, Schwartz (1997) and Manoliu and Tompaidis (2002) employ a state space 
representation to estimate a number of mean-reverting models with the Kalman 
filter, Brigo et al. (2007) apply MLE for the Merton and VG model, Seneta (2004) 
uses moment matching to fit the VG process, and Brooks and Prokopczuk (2013) 
employ the Markov chain Monte Carlo method to fit various stochastic volatil- 
ity models. While, I am sure, the list can be continued for a while, none of the 
above methods appears satisfactory for our purposes as they are either not general 
enough to work for all our processes or require mathematically involved adapta- 
tions in each case. 

Instead, I propose to minimise an appropriate measure of distance between the ker- 
nel density estimate of historical returns and the analytical probability density of 
the model under appropriate constraints and with a suitable algorithm. The idea is 
to choose the process parameters in such a way that the model implied return den- 
sity mimics the empirical return distribution as closely as possible. As suggested, 
for instance, by Schoutens (2003) or Brigo et al. (2007), this is what characterises 
a good model fit. For easier reference, I will call the proposed method Algorithmic 
Probability Density Fitting (APDF) in the remainder of this document. To famil- 
iarise the reader with the idea of APDF, we will next consider the choice of an 
objective function and an appropriate algorithm. Thereafter, the calibration results 
are benchmarked against MFE estimates and some implementation details are con- 
sidered. 

Plausible objective functions include, for instance, the root mean squared error, 
the sum of squared errors, the average percentage error, and the sum of absolute 
errors. Since, however, a measure involving squaring the distance penalises large 
deviations of the model from the empirical distribution more than many small ones 
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(with equivalent sum of absolute errors), the sum of squared errors fits my desire 
to replicate the distribution shape as accurately as possible. While it could be con- 
sidered to add some weight function to put more emphasis on deviations in the 
distribution tail, I leave such experiments to further study as the choice of a sen- 
sible weight function deserves some separate attention. The quantity we want to 
minimise is given by 

N 

SSE = £ (f Kernel (0 " / Model (»)) , (4-2.1) 

i=l 

where SSE denotes the sum of squared errors, f Kernel is the kernel density estimate 
of past returns for a given commodity calculated according to (3.2.13), and f Model 
is the model implied probability density. Where the model density is not explicitly 
available, it is derived from the characteristic function of the process. The SSE 
is computed on a grid of N = 2 16 points, where an exponential of two is chosen 
to speed up the Fourier transform algorithm when the density is computed from 
the characteristic function. 44 The grid is chosen such that it encompasses the en- 
tire spectrum of historically observed returns. Note that in subsequent analyses of 
model fit, I will report the Average Squared Error (ASE), which is the SSE divided 
by the number of grid points N. This simply facilitates the comparison of results in 
case a different number of grid points is used in the density calculations. 

Next, to minimise the objective function there is a choice between direct search 
methods (Nelder-Mead), Newton based algorithms (Levenberg-Marquardt, L-BFGS, 
Sequential quadratic programming), and stochastic optimisers (Differential Evo- 
lution (DE), Simulated Annealing (SA)), which are globally convergent, i.e. do 
not get stuck in local minima (Kienitz & Wetterau, 2012). As experiments reveal 
that local optimisers do not converge to a sensible solution and since computation 
time is not a highly critical argument in our context, it follows logically to choose 
among the stochastic optimisers, where the solution is less dependent on the choice 



44 Sec. for instance, Kienitz and Wetterau (2012). 
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of initial parameter values and, thereby, gives us somewhat greater flexibility in the 
application to different models. Further experiments show that DE is even more 
reliable than SA, albeit lower computation speed. Eventually, I decide to use DE 
in the form originally proposed by Stom and Price (1995) using the open source 
DeMat toolbox. 45 To further improve calibration results and save computation 
time, I combine DE with Pattern Search (PS), a direct search Matlab routine for 
constrained optimisation problems. In other words, once DE has finished compu- 
tation, the parameter results are passed on to initialise PS for further optimisation. 
A scheme of this kind is also known as hybrid optimisation, meaning that a global 
optimiser is used first to circumvent local minima in the objective function fol- 
lowed by a faster local optimiser to fine-tune results (Kienitz & Wetterau, 2012). 
Please note that I use PS instead of other, potentially faster hybrid optimisers as 
it can cope with discontinuities in the objective function. While a detailed discus- 
sion of the algorithms is beyond the scope of this paper, let me say that the logic of 
DE relies on the principles of inheritance, mutation, recombination, selection, and 
crossover. Roughly speaking, the algorithm starts out with a an initial population 
of candidates that are evaluated based on their fitness to minimise the objective 
function. This population is subsequently mutated in a way that the fittest candi- 
dates evolve and an optimal value for the objective or fitness function is reached. 46 
In order to verify next that the method yields acceptable calibration results, I com- 
pare parameter estimates for GBM, the Merton, and VG model under MLE and 
APDF. The choice of these processes is motivated by the desire to have a range 
of diverse models, which admit the calibration by a popular alternative calibration 
method such as MLE for benchmarking. Since the technicalities of MLE calibra- 
tion are not of primary interest in our context, details are summarised in appendix 
B. The graphical representation of calibration results shows that both MLE and 

45 The toolbox is available from http: //wwwl . icsi.berkeley.edu/-storn/code.html. 

46 In subsequent calculations the initial population is set to 20 or 30. depending on the number of parameters to be 
estimated, the mutation scaling factor is set to 0.8, the crossover probability is 0.85, and the number of generations 
is set between 500 and 1000 depending, again, on the number of parameters to be estimated. PS is run with default 
settings. Note that this set-up is conservative in the sense that also fewer population members and mutation generations 
should lead to an optimum. For details see Storn and Price (1995). 
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GBM Merton model VG model 




Figure 14: Comparing MLE and APDF. The plots show the empirical kernel density estimate (green) 
and the analytical return distribution for each model based on MLE parameter estimates (light grey) 
and APDF (dark grey). Parameter estimates (MLE/APDF) for weekly crude oil log returns over the 
full sample period (02/08/1993 - 30/12/2013) for GBM: /j = 0.0016/0.0056, a = 0.0531/0.0464. Cor- 
responding parameter estimates for the Merton model: /j = 0.0016/0.0016, a = 0.0428/0.0366,^/ = 
—0.0291/ — 0.0565 , Gj = 0.1255/0.01, A, = 0.0598/0.3271. Parameter estimates for the VG model: /j = 
0.016/0.0299,0 = 0.0506/0.0481,0 = -0.0144/ - 0.0290, v = 0.6011/0.2467. The empirical kernel and 
analytical (transition) densities are evaluated based on the above stated parameters on a grid of 1000 points. 



APDF parameter estimates yield a close fit to the historical return distribution, 
however, the fit generated by APDF is evidently even closer, especially for the VG 
process. This insight is quantified in table 6, where we can observe that the ASE 
between historical and model implied densities is lower if parameters are calibrated 
with APDF instead of MLE. 47 From these results, we can conclude that APDF can 
be considered an appropriate technique for empirical parameter estimation in the 
next subsection. 

Before we turn to the goodness of fit results, however, let me outline three minor 
implementation details necessary to assure that the APDF delivers sensible param- 
eter estimates. First, in the case of the Fleston and Bates model, I constrain the 
long-term variance 0 to equal the initial variance vo in order to keep the return 
density independent of time. If, for instance, the variance process was allowed to 
decrease to a lower long-term level, the return distribution would become more 
peaked as time progresses. Since this time dependence cannot be captured by the 
objective function, calibration results would be potentially unreasonable. Note that 



47 Note that this conclusion does not change if other distance measures such as the sum of absolute errors or average 
absolute errors are used, which were not minimised in the APDF objective function. 
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Table 6: 

MLE and APDF model fit 

For each calibration method and model, the parameter sets described in fig. 14 are used. In the Analytical 
columns, the ASE is calculated between the empirical kernel density estimate and the model density given 
by its probability density function as in fig. 14. For the Kernel columns, 100,000 returns are generated for 
each model and their kernel density estimate is compared to the empirical one. This comparison is added to 
double-check the reliability of parameter estimates in Monte Carlo applications. 





GBM 




Merton 






VG 




Analytical 


Kernel 


Analytical 


Kernel 


PDF 


Kernel 


MLE 


0.142 


0.168 


0.096 


0.068 


0.076 


0.066 


APDF 


0.038 


0.039 


0.013 


0.014 


0.007 


0.007 



the necessity of this constraint is a weakness of APDF as a calibration methodol- 
ogy for models exhibiting such time dependence in the return distribution. Second, 
I impose the Feller property on the Heston and Bates model. 48 This ensures that 
the processes, in particular the variance process, cannot hit zero. Clearly, prices of 
zero or deterministic prices without uncertainty are not meaningful projections of 
the future. Third, in the Merton, Bates, and NIG-CIR model, I constrain the drift 
to equal the historic mean return. For the former two models this is helpful as the 
algorithm could otherwise compensate an overestimated mean with overestimated 
negative jumps. Controlling the location of the distribution eliminates this poten- 
tial problem and improves convergence, particularly in the NIG-CIR model. While 
this kind of instability issue can always occur if two or more parameters control 
the distribution shape in a similar way, this is an unavoidable evil in the other cases 
and needs to be controlled via plausibility checks. 

4.2.2 Goodness of fit 

In this subsection, it is discussed which model is likely to yield the best fit to his- 
torical commodity returns. To do so, we will draw on the previously developed 
logic and compare the ASE between the empirical distribution estimate and the 
model implied distribution with APDF-fitted parameters. Since complex models 
with many parameters such as the Bates or NIG-CIR process generally outper- 

48 According to Kienitz and Wetterau (2012) the condition holds if 2 k0 > V 2 . 
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form simpler models in in-sample tests, but not in out-of-sample tests (Bakshi, 
Cao, & Zhiwu, 1997), two different test set-ups are created. First, an in-sample 
test, where the ASE is calculated with historical density and parameter estimates 
based on the full sample period from 1993 to 2013. Second, an out-of-sample test, 
where model parameters are calibrated over the period 1993 - 2004 and the result- 
ing model distribution is compared to the empirical density of the period 1995 - 
2013. Clearly, the choice of calibration and analysis period in the out-of-sample 
test is drastic as several extreme events occurred during the analysis period (ascent 
of prices after 2005 and financial crisis), which are very different from events ob- 
served during the calibration period. An additional way of out-of-sample testing 
could involve a leave-one-out cross validation as suggested by De Laurentis et al. 
(2010), however, this exercise is left to further study as the main purpose is equally 
well accomplished by the here adopted simple alternative. This is to understand if 
the nuances in the shape of return distributions that more complex models are able 
to capture persist over time so that it is worthwhile to account for them. In order to 
conserve space, I report the calibrated process parameters after some sanity checks 
in appendix A4. 

The results, shown in table 7, can be analysed in a number of ways. Considering 
first the in-sample results (setting one), we may notice that GBM leads to the worst 
fit in all cases and the Bates model to the best in most cases. Only for nickel, zinc, 
and silver the VG process fits better to past data and in the case of corn, the NIG 
model ranks best. When judged by the mean ASE (in brackets) across commodi- 
ties, we have the following rank order from best fit to worst: Bates (0.0059), VG 
(0.0104), NIG (0.0108), NIG-CIR (0.0156), Merton (0.0260), Heston (0.0270), 
GBM (0.0962). As the Bates model is with nine parameters in some sense the 
most flexible one, the results are generally little surprising. However, I would have 
expected the NIG-CIR process to have a better fit than the stand-alone NIG pro- 
cess. However, while this expectation is not confirmed on average, it is fulfilled in 
the cases of oil and cocoa. When analysing the results by commodity, it is evident 
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Table 7: 



Average squared fitting errors 

For each resource and stochastic process, the ASE is reported for two different settings. Setting one 
refers to the in-sample test, where parameters are calibrated and model fit is analysed over the same pe- 
riod: 02/08/1993 - 30/12/2013. Setting two is the out-of-sample test with calibration period 02/08/1993 - 
27/12/2004 and analysis of fit during 03/01/2005 - 30/12/2013. The parameter estimates corresponding to 
this table are obtained from APDF and are reported in appendix A4 and A5. 



Average squared fitting errors 



Resource Setting GBM Heston 



Oil 


i 


0.038 


0.008 




2 


0.054 


0.025 


Aluminum 


1 


0.143 


0.043 




2 


0.499 


0.382 


Copper 


1 


0.075 


0.029 




2 


0.255 


0.221 


Lead 


1 


0.144 


0.051 




2 


0.756 


0.627 


Nickel 


1 


0.044 


0.006 




2 


0.121 


0.080 


Zinc 


1 


0.093 


0.024 




2 


0.992 


0.892 


Gold 


1 


0.349 


0.087 




2 


2.023 


1.639 


Silver 


1 


0.078 


0.024 




2 


0.366 


0.292 


Wheat 


1 


0.044 


0.027 




2 


0.273 


0.230 


Corn 


1 


0.087 


0.019 




2 


0.269 


0.196 


Soybeans 


1 


0.073 


0.017 




2 


0.189 


0.120 


Sugar 


1 


0.031 


0.006 




2 


0.032 


0.011 


Coffee 


1 


0.020 


0.006 




2 


0.123 


0.099 


Cocoa 


1 


0.127 


0.032 




2 


0.246 


0.133 



Merton 


Bates 


VG 


NIG 


NIG-CIR 


0.013 


0.006 


0.007 


0.009 


0.007 


0.030 


0.025 


0.024 


0.026 


0.024 


0.054 


0.012 


0.021 


0.028 


0.039 


0.453 


0.365 


0.365 


0.373 


0.376 


0.035 


0.008 


0.012 


0.012 


0.028 


0.213 


0.195 


0.202 


0.201 


0.191 


0.063 


0.005 


0.016 


0.020 


0.032 


0.650 


0.570 


0.575 


0.583 


0.590 


0.008 


0.002 


0.002 


0.002 


0.003 


0.084 


0.080 


0.075 


0.075 


0.076 


0.026 


0.010 


0.008 


0.010 


0.016 


0.886 


0.868 


0.856 


0.862 


0.867 


0.066 


0.011 


0.039 


0.024 


0.033 


1.534 


1.485 


1.569 


1.544 


1.522 


0.024 


0.005 


0.005 


0.005 


0.012 


0.280 


0.259 


0.264 


0.263 


0.269 


0.008 


0.005 


0.008 


0.006 


0.008 


0.229 


0.227 


0.228 


0.225 


0.230 


0.013 


0.004 


0.003 


0.003 


0.007 


0.174 


0.173 


0.176 


0.175 


0.181 


0.017 


0.006 


0.007 


0.006 


0.013 


0.106 


0.104 


0.110 


0.112 


0.114 


0.006 


0.004 


0.005 


0.004 


0.006 


0.022 


0.008 


0.009 


0.008 


0.012 


0.006 


0.001 


0.003 


0.003 


0.006 


0.101 


0.095 


0.099 


0.098 


0.099 


0.026 


0.007 


0.010 


0.020 


0.010 


0.115 


0.108 


0.114 


0.131 


0.109 
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that gold has the highest mean ASE (0.087) and coffee the lowest (0.006), indicat- 
ing that processes have on average the greatest difficulties to mimic the shape of 
the gold return distribution while it works best for coffee. 

Turning to the out-of-sample setting, we may notice that the fit generally wors- 
ens. While the mean model ASE per mean commodity is 0.027 in setting one, it 
rises to 0.353 in setting two. The extent to which model fit worsens depends, how- 
ever, on the commodity we consider. For instance, the mean ASE for sugar rises 
only relatively little from 0.009 to 0.015, which is in sharp contrast to the case 
of gold, where the mean ASE jumps from 0.087 to 1.617. This can be explained 
by the changing gold price dynamics, which exhibit an almost steadily increasing 
level of volatility over time (see appendix A2). Although, this was not enough 
to reject stationarity of gold returns previously (see table 4), even this change 
of return distribution between calibration and analysis period suffices to worsen 
model fit considerably. With respect to the rank order of models, it is somewhat 
surprising and almost in contrast to the argument that processes with many pa- 
rameters are easily overfitted 49 that the Bates model retains the best average fit 
also in the out-of-sample test. However, the relative advantage over the second 
best process shrinks. 50 Also, the Bates model delivers only the best model fit for 
50 per cent of the commodities compared to 71 per cent in setting one. For the 
purpose of completeness, we have the following rank order in setting two: Bates 
(0.3258), NIG-CIR (0.3329), VG (0.3334), NIG (0.3340), Merton (0.3483), He- 
ston (0.3534), GBM (0.4428). Please note that these figures are, again, average 
ASEs over all commodities for each model. 

In summary, the analysis shows that in both test settings advanced stochastic pro- 
cesses with more parameters relative to GBM lead to considerably lower fitting 
errors. While it is not uniformly clear across commodities which model is the 
most favourable one, particularly the Bates, VG, and NIG process appear gener- 



49 See, for instance, Bakshi et al. (1997). 

50 In setting one, the mean ASE is 42.6% lower for the Bates model than the second best VG process. This advan- 
tages shrinks to 2. 1 % vis-a-vis the second best NIG-CIR process in setting two, the out-of-sample test. 
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ally good choices. Accounting furthermore for potential difficulties in calibrating 
the high number of parameters in the Bates model, the VG or NIG model resemble 
almost equally powerful and easy to handle alternatives. 
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5 Capital budgeting implications 



Having identified considerable differences in the ability of stochastic processes 
to replicate the properties of historical price series, we do not yet know to what 
extent process choice also influences the valuation of capital investments. While 
Lo and Wang (1995) show that the choice between GBM and a mean-reverting 
Ornstein-Uhlenbeck process can affect financial call option values in the order of 
5%, Tsekrekos et al. (2012) identify much larger variations in a more complex 
real option investment ranging between -40 % and 25%. However, their analysis 
is entirely based on comparisons between the classical mean-reverting commodity 
price processes of Schwartz (1997) so that an analysis of valuation implications 
based on the different perspective advocated in this paper is yet unavailable in aca- 
demic literature. As a result, two questions need to be answered in this section. 
First, does the choice among previously discussed stochastic processes influence 
the valuation result of commodity-linked contingent claims? Second, if valuation 
differences arise, are these equally pronounced across investments with different 
degrees of flexibility or maturities? The latter of the two questions follows from 
the diverse findings of Lo and Wang (1995) and Tsekrekos et al. (2012), which im- 
plicitly suggest that stochastic process choice matters more in complex real option 
applications than in the context of financial options. 

5.1 A stylised investment project 

In order to answer the previously formulated questions, a stylised capital invest- 
ment project with different levels of flexibility is considered. The set-up of this 
project is motivated by the desire to expose the influence of stochastic process 
choice on capital budgeting decisions from a range of different angles. In this 
highly stylised example, a natural resource company has the opportunity to bid on 
the extraction rights for ore bauxite, the basis for aluminium production. Under 
the agreement, the firm is entitled to start production whenever desired and it is 
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estimated that once started, it will take five years to exploit the deposit. The re- 
quired production site at the mine can be set up with negligible time delay after the 
investment decision has been made. Now, to turn this setting into a highly stylised 
valuation case study, the project value will be given as a combination of a static 
NPV, an option to defer investment for five years (decision period), and the option 
to abandon the site for a salvage value (net of administrative accompaniment) once 
production has been started (operating period). For simplicity, it is assumed that 
the operating period is limited to five years. Alternatively, one could fix the total 
quantity of ore bauxite in the mine and let the project last until the resource is fully 
depleted (Brennan & Schwartz, 1985). All inputs and corresponding notation re- 
lated to the investment project are summarised in appendix A6. 

In order to flexibly implement the valuation for different stochastic processes, I 
use the Least Squares Monte Carlo (LSM) technique proposed by Longstaff and 
Schwartz (2001) to estimate continuation values and optimal stopping rules in the 
valuation. This is in line with Gamba (2003) and Tsekrekos et al. (2012), who 
show that LSM can be successfully applied not only to American financial options, 
but also to potentially more complex capital investment projects. On this basis, I 
compute three quantities. First, the static project NPV based on the risk-neutral 
expectation of future aluminium prices if the project is started immediately. Sec- 
ond, the NPV with the flexibility to abandon the project for a salvage value within 
the five operating years but, again, assuming an immediate project start. Third, the 
option value to defer investment for up to five years. From these quantities, the 
separate option values and the total project value including all flexibility can be 
calculated. 

To begin with, the project’s operating cash flow at time t based on stochastic pro- 
cess m and commodity price path CO is given by 

CFo (t , S m (r , CO) ) = qr ■ (S m (t,(S)) - c)- K - Tax (5.1.1) 

Tax = [x (qr ■ (S m (t , ©) - c) - K) , 0] + . 
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Please note that superscript + denotes the maximum of a set of values. The cash 
flow calculation is similar to the one used in the classical natural resource invest- 
ment example by Brennan and Schwartz (1985), however, I omit several details 
such as the adjustment of costs to inflation to avoid any unnecessary complexity 
that distracts from the main idea of this section. Next, I calculate the gross project 
value (gross of the initial investment) for each possible starting date with and with- 
out incorporating the possibility to abandon the project for a salvage value (T|). As 
future commodity prices are uncertain and not part of the investor’s information 
set, the LSM method is used in the backward recursion to estimate the expected 
project value in the subsequent period. Beginning in the last of the five operating 
years Tj of the project when started in year j, where 0 < j < 5, the project value in 
Tj is given by 



i.e. the maximum of the operating cash flow in the last period and the salvage 
value received upon abandonment. 5 1 Now, for t / = T: — 1 , . . . , 0, the project value 
in time t is given recursively by 



where <f> denotes the expected project value in the subsequent period under the 
risk-neutral, equivalent probability measure Q. More formally, this quantity can 
be written as 



In other words, the project value in each period is either the sum of the operating 
cash flow received in that period and the sum of expected discounted future cash 



PV(Tj,S m (Tj, at)) = [CF 0 (Tj,S m (Tj, co)),r|] + , (5.1.2) 





.ti—tj-'r 1 



51 Alternatively, one could have modelled the salvage value in direct proportion to the remaining quantity of ore 
bauxite in the mine, however, this introduces some unnecessary complexity to the problem. 
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flows until maturity if operation is optimal. Or, if operation is not optimal, the 
project value in period t is simply the salvage value as no future cash inflows 
can be expected once the project is abandoned. To account for the possibility 
of different starting dates due to the option to wait, subscript j was introduced 
to indicate the respective starting time of the rolling five-year operating horizon. 
More specifically, if the project is started in year three, the corresponding present 
value in period three is calculated with j = 3 and Tj = 8, following from the five- 
year operating horizon. However, as stated previously, future aluminium prices are 
uncertain so that (5.1.4) cannot be readily calculated at a given time step (Gamba, 
2003). Instead, the continuation value <3> is embedded in the project value step-by- 
step during the backward recursion according to 

PV (i tj,S m (tj, to)) = {CF 0 ( tj,S m (tj, co)) (5.1.5) 

+ e~ rAl E^ [PV (tj + l,S m (tj + l,<a))],i\} + . 

This is often referred to as the Bellman equation (Tsekrekos et ah, 2012), which 
does not involve all future operating dates tj to calculate <f> as would (5.1.4), but 
only relies on the project value one step ahead. This is estimated according to 
the LSM method by first regressing the project’s discounted value of the subse- 
quent period onto an appropriate basis {L n (t,CF a (t,S))} and subsequently using 
the obtained regression coefficients to project the expected continuation value as 
a function of current operating cash flows evaluated through the set of basis func- 
tions. Omitting some of the previous notation for convenience, we can express <E> 
as 

N 

e~ rAt E f Q [PV (t + 1) | CF 0 (r)] = £>„ (t)L n (CF 0 (t,S))+ e, (5.1.6) 

n = 1 

where <j) is the set of regression coefficients corresponding to the value of each 
basis function through which the project cash flow at a given time and commodity 
path is evaluated. The set of coefficients is chosen such that the residual term 
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e is minimised in a least squares sense (Kienitz & Wetterau, 2012). Given the 
coefficients, the continuation value of a project at time t is given in sparse notation 
by 



&(t,CF 0 (t,S))=§ 1 L 1 {t,CF 0 (t,S))+,...,+§NL N (.t,CF 0 (t,S)), (5.1.7) 

where the first six terms of the family of weighted Laguerre polynomials, denoted 
by L, are used as basis functions (see appendix B5). 52 

Once the backward recursion is completed for each starting date in the decision 
horizon for the case with and without abandonment flexibility ((5.1.2) and (5.1.5) 
are modified accordingly), I obtain two grids of gross project values with dimen- 
sions equal to the number of paths and time steps in the decision horizon. Taking 
the average of project values across all paths in t = 0 and subtracting the initial 
investment, I obtain the corresponding NPV with and without abandonment flexi- 
bility when the project is started immediately. As a last step, the additional value 
arising from the possibility to defer investment can be calculated as a call option 
on the gross project value with maturity equal to the length of the decision horizon 
and a strike price equal to the investment cost adjusted for the time value of money 
(Trigeorgis, 1996). Thus, the value of the option to defer (for each commodity 
price path to) is given by 

V(T,co) = [PV(T,ta)-Ie rT ,0] + (5.1.8) 

V(f, 00) = [PV(f,C0)-/e'\V(f+ 1,C0)] + , t = r— 1,...,0. 

As outlined above, the continuation value, denoted by V (t + 1 , co), is estimated re- 
cursively with the LSM method. 53 The option value to wait is finally computed as 
the average across all paths V in time zero. 



52 See Abramowitz and Stegun (1964) for further reference. Please note also that the choice of other eligible or- 
thogonal basis functions such as Hermite or Chebyshev polynomials converge to similar results if enough terms are 
included (Tsekrekos et al., 2012). As a result, no further experiments are conducted in this direction. 

53 Note that contrary to Gamba (2003) or Tsekrekos et al. (2012), I use only so-called in-the-money paths in the 
estimation of regression coefficients as this improves convergence of results (Longstaff & Schwartz, 2001). 



71 



Before turning to the valuation results, let me kindly refer you also to the com- 
plementary Excel spreadsheet where a fully worked out example is provided for a 
more lively illustration of this investment project. 

5.2 Results 

The results presented in table 8 indicate that the choice of stochastic processes 
has a considerable effect on project valuation. It stands out that the use of GBM 
leads to the lowest valuation estimates throughout. On average, GBM underval- 
ues the project by 14.3% relative to the empirically more appropriate alternative 
processes. One explanation for this finding may result from the classical idea that 
higher uncertainty creates value when payoff structures are asymmetric. In other 
words, when volatility is high, an option holder benefits from more frequent ex- 
treme movements as he can gain when prices evolve in his favour, but has a limited 
downside when the opposite is true (Hull, 2009). As all processes except GBM can 
account for the empirically observed fat-tails, i.e. a higher probability of extreme 
events than presumed under GBM, it is plausible that project values should be 
higher using such processes. This logic is consistent with the observation that val- 
uations are highest under the Bates and NIG-CIR process, which imply the com- 
bined uncertainty arising from jumps and stochastic volatility. 

In order to ascertain that the previously presented results cannot be explained by 
sampling error, I use the same pseudo random number stream in each valuation 54 
and provide upper and lower bounds for the valuation estimates at the 95% confi- 
dence level in panel A of fig. 15. 

Besides the impact of different processes on valuation results, it is of practical rel- 
evance to investigate how optimal investment behaviour depends on the stochastic 
price model. To this end, I depict the critical aluminium price as a function of time, 
above which management should proceed with investment. 55 As we would expect, 

54 Please note that valuation results do not change subject to the choice of a different random number stream. 

55 In mathematical finance, this critical price is often referred to as the free boundary that separates early exercise 
from continuation region in American-style derivatives (Battauz & Ortu, 2009). 
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Table 8: 



Capital investment project: Valuation results 

The parameters underlying the valuation results are reported in appendix A4 and A6. NPVq denotes the static 
project net present value without any flexibility. 02A refers to the option value from abandonment. NPV \ 
resembles the project NPV including this abandonment flexibility. 02 W denotes the option value to wait. 
NPV 2 is the project value including both the flexibility to abandon and to wait. CRR refers to the valuation 
results obtained from the binomial lattice of Cox et al. (1979) with parametrisation corresponding to GBM 
in the Monte Carlo simulation. The CRR valuation is added as a reference point and to confirm the reliability 
of the LSM algorithm in this setting. Please note that all valuations are performed using the same global 
pseudo random number stream for better comparability. Results are reported in US$ ’000s. 





CCR 


GBM 


Heston 


Merton 


Bates 


VG 


NIG 


NIG-CIR 


NPV 0 


-256.92 


-256.95 


-266.37 


-266.08 


-266.64 


-268.33 


-269.37 


-268.30 


+ 02A 


294.60 


293.40 


316.72 


317.59 


326.23 


320.55 


324.11 


324.37 


= NPVi 


37.68 


36.45 


50.35 


51.50 


59.58 


52.23 


54.74 


56.07 


+ 02W 


365.29 


368.55 


398.65 


400.98 


428.15 


403.96 


408.85 


412.70 


= npv 2 


402.96 


405.00 


449.00 


452.48 


487.73 


456.19 


463.60 


468.77 



this threshold is declining in time for all processes. Since the decision to invest at 
each step resembles a comparison between the expected value of waiting and im- 
mediate investment and since the value of waiting declines as we move closer to 
the end of the decision horizon, already a lower commodity price (and correspond- 
ing operating cash flow) is sufficient for immediate investment to be optimal. In 
our context, one may observe in panel B of fig. 15 that critical prices are relatively 
aligned, but slightly higher for GBM than the remaining models. This suggests 
that GBM may not only undervalue the project but also lead to overly conservative 
investment decisions relative to the value maximising policy implied by empiri- 
cally more appropriate processes. 

The flipside of this consideration is the probability of investment as a function of 
time. Similarly to Schwartz (2004), I compute this quantity as the percentage of 
paths, for which investment is optimal at each point in time (without conditioning 
on prior exercise). In line with previous arguments, it is evident from panel C of 
fig. 15 that under GBM, management would often defer the project start more than 
optimal. 

While the previous paragraphs have already confirmed the importance of stochastic 
process choice, it is yet unclear whether this conclusion is equally valid for diverse 
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A: Valuation overview 
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450 
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B: Critical aluminium price 



0.06 



0.04 



0.02 



i * 



■ 



G H M B V N NC 
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D: Project complexity 
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Figure 15: Capital investment project: Valuation analysis. G: GBM, H: Heston, M: Merton, B: Bates, V: 
VG, N: NIG, NC: NIG-CIR. Panel A depicts the project valuation in US$ ’000s including the flexibility to 
defer and to abandon. Also shown are valuation confidence bounds (95% level). Panel B shows the critical 
aluminium price in US$ ’000s, above which it is optimal to invest and below to wait. This is plotted for 
the 5 year decision horizon. Panel C illustrates the fraction of paths for which investment is optimal. In the 
last step of the decision horizon (not shown for more convenient viewing), investment is optimal in 55% of 
cases with little variation across processes. Panel D depicts the percentage deviation of the value of different 
project components for all processes relative to GBM. Here, NPVq is the static NPV, 02A refers to the option 
to abandon, 02AW is the value of the option to wait on the project, which can be abandoned, and NPV \ is 
the project NPV including all flexibility. 



projects with different degrees of complexity or flexibility. To answer this second 
question in more detail, panel D of fig. 15 provides the percentage difference of 
project value for all stochastic processes relative to GBM when different degrees 
of project flexibility are considered. The idea is that if complexity exacerbates 
the impact of process choice on valuation results, percentage differences in project 
value should gradually increase for different processes when more optionality is 
present in the project. Indeed, it is evident that valuation estimates are relatively 
similar when only the uncertain, static NPV is estimated, however, as more op- 
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tionality is added, estimates diverge. This suggests that an appropriate commodity 
price model is essential for valuation accuracy in complex investment decisions, 
however, it may even be negligible in simple cases. To confirm this logic, I also 
compute the price of a simple at-the-money call option on one ton of aluminium 
with corresponding five-year maturity and observe that all estimated prices differ 
only in the order of three to eight per cent, i.e. much less than in the previous case 
where payoffs were more complex. 

To sum up, these findings are generally in line with earlier studies of Lo and Wang 
(1995), Casassus and Collin-Dufresne (2005), Tsekrekos et al. (2012), and Brooks 
and Prokopczuk (2013), who similarly indicate the relevance of stochastic process 
choice in the valuation of financial and real assets. While the magnitude of valua- 
tion differences ranges between -40% and 30% across studies, it has become clear 
that this can be explained by diverse levels of uncertainty and flexibility, as well as 
the length of investment horizons and the underlying commodity. 
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6 Conclusion 



This thesis aimed at providing practically relevant recommendations for the choice 
of stochastic commodity price models in capital investment valuation. The anal- 
ysis of empirical price dynamics across a basket of 14 commodities led to the 
conclusion that popular mean-reverting processes and GBM are in several respects 
inconsistent with empirical data. In turn, the merits of alternative models known 
from financial options pricing were analysed in the application to commodity mar- 
kets. In this context, a practical and flexible calibration scheme was proposed that 
allowed me to assess the relative goodness of fit of each process to a range of 
commodities. It was found that all new models deliver a considerably better fit 
to historical commodity price data than GBM. Among the tested models, the best 
fit is on average achieved by the Bates and VG process, however, accounting for 
the unfavourably high number of parameters in the Bates model, the VG process 
resembles a sensible compromise between accuracy and complexity. To gauge the 
extent to which the choice among different stochastic processes deserves careful 
attention in practice, a stylised investment project was introduced. The analysis 
revealed that in this example stochastic process choice may affect project values in 
the order of 20%. Moreover, it was found that GBM may not only underestimate 
project value, but also lead to excessive investment deferral under the optimal deci- 
sion rule relative to empirically more appropriate processes such as the VG model. 
While the valuation example was exclusively based on aluminium price dynamics, 
it can be expected that empirically more skewed or fat-tailed commodities such as 
soybeans or coffee lead to even higher deviations from GBM-based estimates. 
These findings may complement existing literature in three main aspects. First, 
they are based on a relatively broad dataset and bottom-up statistical approach to 
stochastic process selection with practical recommendations. Second, this study 
is, to my knowledge, the first to introduce Levy processes to the world of real 
options. Third, stochastic models are not considered in isolation, but comprehen- 
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sively evaluated in their impact on investment valuation and optimal managerial 
decision making. 

However, several topics follow from these insights for future study. In particu- 
lar, the universe of stochastic processes considered here is only a small subset of 
modelling possibilities and can be regarded merely as an entry point for future dis- 
cussion. For instance, the role of stochastic interest rates outlined by Schulmerich 
(2010) or seasonal price components as mentioned by Brooks and Prokopczuk 
(2013) have not even been touched upon. Furthermore, the considerations of this 
thesis could be extended from a single source of uncertainty to the joint distribu- 
tion of several risk factors that determine investment valuation and optimal deci- 
sion making. With respect to the proposed APDF calibration scheme, it may be 
useful to assess to what extent the more forward looking calibration to financial 
option prices affects parameter estimates and valuation results. Finally, the LSM 
implementation can be refined by computing regression coefficients and project 
values based on different random number streams to eliminate a potentially biased 
valuation estimate (Kienitz & Wetterau, 2012) and antithetic or control variates 
could be applied to reduce Monte Carlo variance in the study of stochastic pro- 
cesses (Battauz & Ortu, 2009; Tsekrekos et ah, 2012). 

While this list can be continued in several directions, the overriding principle in fu- 
ture research should mainly reside in the practical application of real options theory 
to support managerial decision making in an increasingly uncertain environment. 
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Appendix 
Appendix A 




Al Co Le Ni Zi 

Industrial metals 



B: Illiquidity over time 




Starting date 



Al: Liquidity of commodity prices. Al: Aluminium, Co: Copper, Le: Lead, NI: Nickel, Zi: Zinc. Panel 
A shows the Zeros metric for daily spot and 3-months future prices from the LME for the period 02/08/1993 
- 30/12/2013. Panel B illustrates the development of Zeros for daily spot and 3-months future prices for 
aluminium only. The x-axis denotes the respective starting date for one-year revolving evaluation periods. 
For convenient viewing only 100 data points are plotted. Note that all non-trading dates have been cleared 
from the data prior to the analysis. 
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A: Energy 




B: Industrial metals 




C: Precious metals 




D: Agriculture 




Time 

A2: Historical price evolution. Al: Aluminium, Cc: Cocoa, Cf: Coffee, Cn: Com, Co: Copper, Go: Gold, 
Le: Lead, Ni: Nickel, Ol: Crude oil, Si: Silver, So: Soybeans, Sp: S&P 500, Su: Sugar, Wt: Wheat, Zi: 
Zinc. All data is presented for the full sample period: 02/08/1993 - 30/12/2013. 
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A5: Out-of-sample APDF: Calibrated parameters (1) 

Al: Aluminium, Cc: Cocoa, Cf: Coffee, Cn: Corn, Co: Copper, Go: Gold, Le: Lead, Ni: Nickel, 01: Crude oil, Si: Silver, So: Soybeans, Su: Sugar, Wt: 
Wheat, Zi: Zinc. Calibration period: 02/08/1993 - 27/12/2004. Data frequency: Weekly. 
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0.0055 0.0076 0.0539 0.4613 



A5: Out-of-sample APDF: Calibrated parameters (2) 

Al: Aluminium, Cc: Cocoa, Cf: Coffee, Cn: Corn, Co: Copper, Go: Gold, Le: Lead, Ni: Nickel, Ol: Crude oil, Si: Silver, So: Soybeans, Su: Sugar, Wt: 
Wheat, Zi: Zinc. Calibration period: 02/08/1993 - 27/12/2004. Data frequency: Weekly. 



o cn o in 
o o o in 
© © © © 



CN CO NO 
_ co co 
CN NO CN — 
© © © co 
dodo 



o o o co 

dodo 



o CO o CO 
© © © oT 

dodo 



O co O C"- 
© © © or 
dodo 



O O O CN 

dodo 



O cO O or 
© © © oT 

dodo 



o o o oT 

dodo 



o cn o in 

© © © oT 

dodo 



o or o on 
o o o co 
dodo 



o o o or 

dodo 



ON CO h i- 

O ON. — o 
O CN O 1/1 



o o o o 



O O O CN 

dodo 



CO OO O 00 
NO On in 
oT CO t"- 
© © — 
dodo 



3. 0 O > 

0A 



O O NO o 



co O t"- o 
o -h or no 

CN in CN ON 
O ON CO p 

o »n no d 



co on o in 
oo t — r~~~ oo 
O 00 CO NO 



o or or o 



oo o cn oo 
o r- o oo 

O CO t""- CO 

o co oo o 



cn t — on r— 

CN t — t — CO 



O NO CN O 



or in on cn 
or co co in 

P O 00 NO 
P J-; p p 
d co of d 



oooNoor"- 

h - O OO o 

— i r- in no or 

t"~- O O NO 



no r- oo co © © t-- 

o on o in t — O O 

o or on in < o no 

o ^ h cn oo o co 

o no cn o o d d 



— i on in oo o in o 

c - or in - o cn 

O VO On h is in Oi 

O of ■'t 't 'O O' 

o h in o h of o 



no oo no or in 

NO OO ON CO or 

O CO no CN On 

co CN or on On 



co no o in 

O OO CO CO 

o on or ' 



in r~~ 

CO NO 
NO OO 



h CO I — t — O 



O ON o o 



i^incnovo-- 
ON 00 -- -H (s ON 
on o — or n co 

• oo ^ ON oo p- 



O CN CN CN 



o in or co o no or 



o cn o o 

d cn cn o d o o d cn cn d 



cn or on or o co no 

O of ■ 00 O O 00 

o or in cn o on co 

o ^ in o o in p 

© no © © © © © 



© t- © in 
O — I CO — 
O or CO CN 



O CN co O 



CO t — CN oo 
O CN OO CN 
^ oo or 
o r- cn o 



no co co On 



O 00 ON O 

o o d 



cn no o in 
cNOinc-'- 
— * in in or 



©OOCNO OONCNOOOoT 



in I 1 co CN MD h 
© o on on no o »n 

© co — oo <n o no 

© ^: © © on © or 

© © © © © © © 



On co >n NO t — © 00 

rH no OO O CO O VO 

© in On in ON © CN 

© ,-j! CN — < 00 © © 

© r- o © © © © 



no © oo © no © no 

^ © © © CO © I"- 



© © © 



00 ON CO CO © 

© oo or or © 



© ON CN © © CN © 



oo no > >n © 

© or NO 00 CN 

© ON CO CN CO 



•n co co no 
r- or ot- no 
or of CN 
p CN < 

o r- O' o 



3. S CO. CO 

DIN 



CN OO -h CO O CN 

co no fs oo o or 

CN CN NO © © 

in d d d d 



3. S co. to << X p- 

HI0-OIN 



86 



A6: Capital investment project: Valuation parameters 

The table reports input parameter values for the stylised capital investment project used in the comparison of 
stochastic processes. The calibrated parameters for each process and resource can be found in appendix A4. 
Since calibrated parameters for the commodity price are in weekly units, the annual inputs shown below are 
converted to weeks where necessary. A simplified counting convention is used (50 weeks/year). 



Project inputs: 


Project start: 


02/01/2014 


Total time horizon (T): 


10 years 


Maximum operating time: 


5 years 


Decision time horizon: 


5 years 


Aluminium production quantity (qr): 


300 tons/year 


Initial aluminium price (SO): 


$1776. 8/ton 


Variable production cost (c) 


$ 1700/ton 


Fixed costs (K): 


$65, 000/year 


Tax rate (x): 


30% 


Risk-free interest rate (r): 


2% 


Net convenience yield of aluminium (q): 


0% 


Salvage value (r|): 


$20,000 


Discretisation steps per time path 


100 


Number of paths 


150,000 
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Appendix B 

Bl: Expected value and variance for GBM and the Vasicek model: 

This part of the appendix shows how process mean and variance are calculated to 
produce fig. 3. The corresponding process discretisation and MLE implementation 
also used in this analysis are shown in the subsequent paragraph. 

Considering first the case of GBM, Brigo et al. (2007) show that expected value 
E [S' (t) \S (0)] and variance V [5 1 (t ) IS 1 (0)] are given by 

E [S 0)] = S (0) exp 

V [5 (t )] = S (0) 2 exp (2 /jt) ^exp ■ 

In the case of the Vasicek model or Ornstein-Uhlenbeck process, expected value 
and variance are, for instance, reported by Brigo et al. (2007) and are given by 

E [5 (t )] = + (S (0) - /j) exp (-fa ) 

V[S(t)] = ^ l ( 1 - ex P(- 2 ^)), 

where /j denotes the long-run mean level of the process and A, the mean reversion 
speed. 

B2: Maximum likelihood calibration: 

The MLE method was used to calibrate parameters in the backtest presented in 
fig. 3 and the discussion on the quality of results obtained from the APDF algo- 
rithm in section 4.2. 

GBM: In the case of GBM, the log-retums are a series of independent and identi- 
cally distributed random variables, each with density /e (r) = / (f;/u;a) (Brigo et 
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al., 2007). The return density is known to be normal and given by 
/(9) = 7V;(^^-^o 2 ^Ar,o 2 A^ . 

The MLE estimate of 0 i.e. the set of parameters is obtained by numerically max- 
imising 

T 

In X(0) = £ log/e (r t ), 

t= 1 

where T denotes the total number of return observations in the sample. The imple- 
mentation is based on Brigo et al. (2007) and available, as all other calculations, in 
the Matlab package provided complementary to this thesis. 

Vasicek: The Vasicek model is calibrated in levels not returns. The log-likelihood 
function for a set of (commodity) prices is given by 

L Gu, K, o) = £ In/ (S (t) S (/ - 1) ;/i; k o) 

i-l 

= — ”ln(2jt) — Tln^) 

- 2^2 L I 5 (0 _ S ( ? “ !) ex P (~^ Af ) ( 1 - exp ( — A. Ar ) ) ] 2 . 

Besides the possibility of numerical maximisation, explicit maximum likelihood 
equations can be derived. 56 Please note also that I verify the equivalence of pa- 
rameter estimates under MLE and OLS regression calibration. 



56 For further reference see http://www.sitmo.com/article/calibrating-the-ornstein-uhlenbeck 
-model/. 
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Merton: The log-likelihood function for the set of historical returns r is given by 



T 



In L(Q) = In (r t ;{i; a; juj;Oj) 



£P(T> = j)M 




which is an infinite mixture of Gaussian random variables, weighted by a Poisson 
probability P (T t = j) = fp (j;XAt) (Brigo et al., 2007). 

Variance Gamma: In the estimation of the VG process parameters, I follow Brigo 
et al. (2007) and use the Matlab function mle to calculate the MLE estimates from 
the model probability density as given by (4.1.18). The numerical optimisation is 
initialised as follows: 



where K1 and SI are the standard measures for excess kurtosis and skewness out- 
lined in section 3.2.1. 

B3: Characteristic function and Fourier transforms 

According to Schoutens (2003), the characteristic function (j) of a distribution X, is 
the Fourier-Steltjes transform of the distribution function F (x) = P(X < x): 




v = K 1 At 

a Slov/At 

0 = — 



3v 




•infty 



(j) Y (u) = E [ exp {iuX)] = 



exp (iux)dF (*) . 
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The characteristic function always exists and, in our context, characterises the re- 
turn distribution of a stochastic process uniquely. In order to calculate the prob- 
ability density on a discrete grid from the continuous characteristic function, the 
fast Fourier transform algorithm ifft in Matlab) is used. From an appropriate grid 
of N points, the characteristic function values x are transformed to the probability 
density X according to 



and i refers to the imaginary unit in complex numbers. The calculation of probabil- 
ity density on a discrete grid of N returns, requires an appropriate transformation 
of the return grid prior to computing corresponding characteristic function values 
and Fourier transform. The implementation of this procedure is available in the 
complementary Matlab file fft_pdf. 

B4: Quadratic exponential scheme for Heston and Bates process 

The discretisation used in this paper directly follows Kienitz and Wetterau (2012, 
p. 297) and involves the following algorithm: 




where 




1. Compute 




m = © + (V (t) — 0) exp (— kA) 
2 _V(t)v 2 e~ Kd f 



2 
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2. For comparison calculate V|/ 





(1 — exp(— kA)) 2 



(0 + (V (t) — 0)exp (— kA)) 



2 



3. Generate u 11 ( 0 , 1 ) 

4. Compute V (t + A) by switching rule 

if \\i < V|/ c . then 

Use a non-central y} distribution for approximation and mo- 
ment matching. 

A: b 2 = 2\|/ 1 - 1 + ^/rVsrl-l > 0 and a = jfp 
B: Let z\ denote a variate of Z\ tv smi) 

C: Set V (t + A) =a (b + zi ) 2 

if \|/ > \|/ c - then 

Use an exponential distribution to approximate the true distri- 
bution. To generate variates from the exponential distribution 
apply the following rule: 



5. Denoting by zi a realisation of Z ~ (0, 1 ), set 

X(f + A) = X (t)+Ko+KiV (t) +K 2 V (t + A) + a/ K 3 V (t) +K 4 V (t + A)z 2 . 
For a grid with equidistant time steps, as used in this paper, Kq , . . . , K 4 are constants 




C: Set V (t + A) = t|/ l(u;p, (3) with 



¥ *(«)=¥ 1 (3) = < 



0 if 0 < u < p 
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and can be calculated according to: 



Kq = 



OK0 , 
- A 



K\ =YiA 




K 2 = Y 2 A 




K3 =YiA(l-p 2 ) 
^4=Y2A(^l-p 2 ] . 



p 

V 

p 

V 



The martingale correction for risk-neutral simulation in Monte Carlo-based valua- 
tion using the LSM algorithm, can be accomplished by replacing Kq with 



*o = - log (M) - {Ki + 0.5K 3 ) V„, 



where the constant M can be computed explicitly when V|/ < y r according to 



M = 



exp 



( Ab 2 a \ 
\—2Aa ) 



\/l — 2Afif 






When vj/ > tj/j, AY is given by 



,<1 



where 



A = K 2 4- O. 5 A 4 = ^ (1 +KY 2 A) -0.5Y2Ap 2 . 

For further reference on the QE scheme, please refer to Andersen (2008) or Kienitz 
and Wetterau (2012) and for the Matlab implementation consider pathJieston for 
the discretisation under the physical probability measure P and patliJiestotun for 
the risk-neutral counterpart. 
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The corresponding implementation of the scheme adjusted for jumps as it is used 
for the Bates model simulation is provided in Matlab file pathLbates. 



B5: Basis functions and confidence intervals in LSM 

Basis functions: The weighted Laguerre polynomials used in this paper were sug- 
gested by Longstaff and Schwartz (2001) and are given by 



The interested reader is also referred to Abramowitz and Stegun (1964) for a com- 
prehensive overview of alternative orthogonal polynomials not considered in this 
thesis. 

LSM confidence interval: These are computed according to Glasserman (2003) 
and Battauz and Ortu (2009). In particular, the radius of the confidence interval 
can be calculated as 



where d is the sample standard deviation of Monte Carlo estimates across paths, n 
is the number of simulated Monte Carlo paths, and 5 is the confidence level (e.g. 
95%). Note that the quantity a n /^/n is often referred to as the standard error of es- 
timates (Glasserman, 2003). More specifically, the volatility of estimates is given 
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by 



<*„ = ^E«-^) 2 , 

where V, is the Monte Carlo estimate of project value for commodity price path i in 
time zero and a is the average value across paths i.e. the Monte Carlo estimate of 
project value. The upper and lower bound for valuation estimates can be computed 
by simply adding or subtracting the radius of the confidence interval to or from a. 
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